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Abstract 


The sum-difference surface current formulation is introduced to treat electromagnetic 
boundary value problems when anisotropic impedances are specified on both sides of a 
surface. It can also be applied to impedance coated bodies. This formulation preserves the 
duality nature of Maxwell equations and carries it over into the algebraic form of the 
integrodifferential operators in the equations for surface currents. Since a 90° rotation is 
equivalent to undergoing a duality transform for an incident plane wave, this particular 
symmetry in the algebraic form of the operators leads to sufficient conditions under which 
the on-axis backscattering of an anisotropic impedance coated scatterer having a 90° 
rotational symmetry is eliminated. The sum-difference formulation is utilized for solving the 
problem of electromagnetic scattering from an anisotropically impedance coated tubular 
cylinder of finite length. The solution has been coded in FORTRAN and tested. Some 


interesting results are presented and discussed. 
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I. INTRODUCTION 


Sometime ago the question was raised: "For electromagnetic boundary value 
problems with specified surface impedances, how can one go from a non-perfectly 
conducting surface on which both the electric and the magnetic equivalent surface currents 
are to be found, to a perfectly conducting surface on which the number of unknowns 1s 
halved [1]?" The answer to this question turns out to be one of algebra. It is well known that 
the impedance specified on the surface of a body separates its interior completely from its 
exterior. Therefore an impedance coated body can always be considered as a hollow volume 
enclosed by an infinitesimally thin shell with surface impedances specified both on the inside 
and the outside of the shell. The inside and the outside of the body can be considered as 
constituted of the same medium and the impressed electromagnetic excitation can be treated 


as continuous across the shell. On the outside surface, there are the equivalent total electric 


current K” and total magnetic current L ’; on the inside surface, there are K and L . For 


an exterior problem, only K* and L” need tobe found; for an interior problem, only K and 


—_ 


L are necessary. A single formulation for solving both types of problems would appear to 


require finding all inside and outside currents therefore doubling the amount of work, but it 
turns out not to be the case because some of the currents are linear combinations of others. 
Furthermore, this formulation holds the key to answering the question posed above. 

Since the shell is infinitesimally thin, from Maxwell equations the radiation to the 


outside and to the inside of the shell can both be given in terms of integrodifferential 
operators on the sum currents K = K’ +K and L=L°+L_. Note that the outside currents 


will not contribute to the radiation in the interior while the inside currents will not contribute 
to the radiation to the exterior of the shell. For simplicity in the description, we consider the 
exterior problem of electromagnetic scattering. By definition, the radiation is the difference 


between the total field and the incident field. On the surfaces of the shell, this definition links 
the incident E and H fields and the difference currents K° -K and L°-L tothe sum 


currents. It is therefore natural to treat the difference currents and the sum currents as the four 


unknowns to be solved instead of the inside and the outside currents. 

The surface impedance on the outside surface of the shell normalized to the medium 
is denoted by Z* and that on the inside surface is Z’. They can be tensors if the impedances 
are anisotropic and may vary from point to point. By forming the sum impedance 
Z=(Z°+Z_)/2 and the difference impedance A =(Z* - Z )/2, the impedance boundary 
conditions provide a set of relations between the difference and the sum currents. It turns out 
that the rank of Z determines how the unknown surface currents are solved. If Z is invertible, 
then the difference currents are linear combinations of the sum currents so that only the 
integrodifferential equations of the sum currents have to be solved. There are only two 
unknowns to be solved for both the exterior and the interior problems under this sum- 


difference current formulation. If Z is rank 0, then Z = 0; the impedance boundary condition 


requires that L be proportional to a 90° rotation of AK. The difference electric current is 


obtained from K after the integrodifferential equation on K is solved. This situation includes 


the case where Z* = Z~ =O when the surface 1s perfectly conducting, thus the result answers 
the question about the transition of the equations from a problem of two variables to one 
which has only a single variable. 

Instead of dealing with an impedance coated body, this thesis presents the sum- 
difference currents formulation of electromagnetic boundary value problems for the 
scattering of an infinitesimally thin surface for which both the inside and the outside currents 
are true unknowns to be found. Extension of this formulation to impedance coated bodies is 
then discussed. 

This formulation preserves the duality nature of Maxwell equations and carries it over 
into an explicit specific algebraic form of the integrodifferential operators in the equations 
for the sum currents. Since, for a plane wave, a 90° rotation is equivalent to undergoing a 
duality transform, this explicit symmetry in the algebraic form of the operators enables us to 
deduce sufficient conditions for eliminating the on-axis backscattering of an anisotropic 
impedance coated scatterer having a 90° rotational symmetry. 


The sum-difference currents formulation is utilized for solving the problem of 


electromagnetic scattering from an anisotropically impedance coated tubular cylinder of 
finite length. The solution has been coded in FORTRAN and tested. Some interesting results 


are presented and discussed. 
In this thesis, the time dependence e ~‘®! is used. E represents the electric field 
intensity divided by the intrinsic impedance of the medium, /p/e; therefore it has the same 


unit as H in amperes per meter. So are the electric and magnetic equivalent surface currents K 


and L. 





II. THE SUM-DIFFERENCE SURFACE CURRENT FORMULATION OF 
ELECTROMAGNETIC BOUNDARY-VALUE PROBLEMS 


A. STRATTON-CHU FIELD FORMULATION AND RADIATION 


On an orientable, piecewise regular surface, whether open or closed, having a surface 
electric current density K and a surface magnetic current density L, we define the Stratton- 


Chu E-field formula as: 


— — — ; 2 — ’ — 
E,_@S,R,D) = = [ RG)GP-7)da,-V [ RG) +V,GF-F,)da, 
Awad S 4% § om 
k ( ia ) 


ay x Lr Gir da 
re [ Led GC da, 


ik|7-7,| 
— a e e ° — —> é 
where 7 is a point which is not on S, k = @,/u €, and G(r-r,) = kiFoF then the 
ear. 


| 


Stratton-Chu H-field formula can be defined as: 
Ce eee) . S. Laeck) 
Lo ik? 
= —— Vx [ K (7, )G(r-r da, + Ae if BMG Gia) aa, (2-2) 
my, LF) °V, G(r-r_)da, 


Note that if S is a closed surface and K and L are the actual total equivalent surface currents 


on S$, then Ee and H.. are respectively the E and H fields at 7 due to all sources inside 
S if 7 is located outside S and vice versa. This is a direct consequence of Maxwell’s 


equations [2] and under this circumstance, the Stratton-Chu formulae are equivalent to 
Maxwell’s equations. On the other hand, unlike the Poynting theorem, the Stratton-Chu field 


formulae over an open surface S are but integrodifferential operators on the tangential vector 


fields K and L over S without any special physical meaning attached. 
To introduce equivalent surface currents on S, the direction of the unit normal vector 7 


on the surface has to be chosen. Adopting the terminology for a closed surface, we can assign 


one side of any orientable surface S as the “outside surface" S$”, albeit somewhat arbitrarily 
if Sis not closed. The outward normal 7 ”* is the unit normal pointing out of this side of § 
and for simplicity, we call 7° the normal of S and denote it by 7. The other side of the 
surface S is the "inside surface" S°. At any point 7 on S, the inward normal nv is -7. Asa 
convention, the fields and surface currents on §* and S” always carry the corresponding 


superscripts (Figure 2-1). 





Figure 2-1. Outside and Inside Surfaces, Normals and the Equivalent Currents. 


Each of the total surface currents K~ and L* on S* consists of two parts: the incident 


current (with the additional superscript "inc") and the scattering current (with the additional 


superscript "sc") corresponding to the incident field and the scattered field on the particular 


side of the surface S: 


ry it a+ y “A apie | Fan ry +, 
K~ =7*x H*> =n*xH™ +4*x HH 
(2-3) 
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hy 
tI 
ty 
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Note that S is infinitesimally thin, hence H"” =H" =H" and BE” =E7'" = E”* on 


— 


Sso that K°"" = -K'" and L'" = -L'"*. Therefore the sum currents K and L on S 


defined below are also the corresponding sums of the scattering currents only: 


2-5 
(. [i720 oe ia 


Since the Stratton-Chu field formulae are linear operators on the surface currents, the 


radiated fields from surface currents on S are determined by the sum currents only: 


E* (nA 7 BGsaK .L ) i Say (Sree a, ) = Ey — 
(2-6) 


He (ea (7 S)K, 0) WES) ere) 


B. CONDITION ON THE CURRENTS IMPOSED BY MAXWELL’S 
EQUATIONS 


As f approaches r* on S*, the tangential components (denoted by the subscript 


"tan") of Eq. (2-6) provide four equations relating the incident fields and the total currents 
on both sides of S through the fact that the incident field is the difference between the total 


and the scattered field: 


EB Sig ERIS) Om 
Eee — BG asa) (2-8) 
Ho = Ho - Ao ea? SBD) (2-9) 
je eles eles (Pag Swhel) (2-10) 


(2-11) 


the difference between Eqs. (2-7) and (2-8) trivially confirms the definition of the sum 
equivalent magnetic current while the difference between Eqs. (2-9) and (2-10) confirms the 
definition of the sum equivalent electric current as both can be deduced directly freien 
Maxwell’s equations. We choose to use the sum of Eqs. (2-7) and (2-8) and that of Eqs. (2-9) 
and (2-10) as the two independent linear combinations out of Eqs. (2-7) through (2-10) to 


link the incident fields to the total surface currents on S* as dictated by Maxwell’s equations: 


il 
=> 
4 

BR 

i 
a 
aa 
= 


where M and N are linear integrodifferential operators on the tangential vector fields K and 
L over S. 
Under any orthonormal coordinates (u, v) over S having i, ¥ as the unit basis vectors 


and with 7 = @ x ¥,a tangential vector field A over S can be written in matrix form as: 


A 


— oy 
Aa 
A 


i | is one of the Pauli spin 


u 


.Thenz x A = = Ou where Oo = 

















Vv u 


matrices. Note that o,* = 1. Using these matrix notations, we can rewrite Eq. (2-12) in the 


following form: 

















0 i0))K°-K| [mM -N|/R] | |Ean 
S oa jem | (2-13) 
-io, O}lFt_f-| |N MUIE At 


ce. IMPEDANCE BOUNDARY CONDITION . 
Maxwell’s equations alone cannot determine the electromagnetic fields completely. 
If S is an open surface, appropriate boundary condition which the fields satisfy on S must be 


specified. It is usually given in terms of the impedance boundary condition, a linear relation 
among the tangential components of the total E and the total fields on S. If S is a closed 


surface, there are two possibilities: One is to specify the electric and magnetic properties of 


the volume within S and require the fields to satisfy regularity conditions within S and be 
linked to the fields outside through boundary conditions across 5S; another is to specify the 
impedance boundary condition on S * for an exterior problem or on’S for an interior 
problem. Note that an impedance boundary condition over a closed surface S completely 
separates the exterior from the interior of S. Therefore, the surface impedance on S” can be 
arbitrary for an exterior problem while that on S * can be arbitrary for an interior problem. 
In this thesis, normalized surface impedances Z* are assumed to be specified on S * whether 
Sis an open or aclosed surface. Note that an open surface S can be considered as bounded 
within the closed surface formed by joining S* and S~. 


The impedance boundary conditions on S* are defined by: 


x n*) = Z* (i* x H*) (2-14) 


or equivalently, in terms of the total surface currents: 


n= x Lo= Z* K~ (2-1 


With the matrix notations for tangential vector fields over S$ in the orthonormal coordinate 
system (u, v) introduced before, we can consider Z* as 2X2 matrices and rewrite Eq. (2-15) 


in the form: 


which can readily be recast into a relation among sum and difference currents: 


K 
(2-175 





= NA ee 
Wh 0 
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with 


Tie Taam (2-18) 


and 


A = Ci a (2-19) 


Eqs. (2-13) and (2-17) are a set of four two-dimensional vector equations to be solved for the 


sum and difference equivalent electric and magnetic surface current densities on S. _ 


D. ALGEBRA OF THE SUM-DIFFERENCE CURRENT EQUATIONS 

The existence and uniqueness of a solution to either the exterior or the interior 
problem specified in terms of the impedance boundary condition have been well established 
[3]. Here we want to investigate how such a solution can be obtained from Eqs. (2-13) and 
(2-17). Clearly Eq. (2-17) defines uniquely the algebraic relationship between the difference 
and the sum currents if Z is invertible. For example, the difference currents can be expressed 


in terms of the sum currents: 














ee a eee 
K' -K 0 id, K 
= = ; R _ (2-20) 
ee = L -t0, 0 L 
where 
fh NTE TS) NTS oe 
jes (2-21) 


“= =I -| 
10,2 A 6,2 Oy 


An equation for the sum currents is obtained by substituting Eq. (2-20) into Eq. (2-13): 




















M -NIIK R ae 
+ R = 
N r r ny ne 
MIL L ae 


ODN 


which can be solved for K and L. Eq. (2-20) in turn enables us to compute the difference 


currents algebraically then split the sum and the difference currents into total currents on S*. 


If Z is not invertible, then the situation is more complicated. Z can either be of rank 


O when Z = 0 or rank | when det Z = 0 but Z # 0. Eqs. (2-13) and (2-17) can be combined 


into an equation for the sum currents K and L: 






































1 rAtoMiye online Za A iver || 5. 
0 -iZo,||N M\\z| |-4 -io,||z| “Jo -iZo,| lq 
am tan 


(293) 


When Z = 0, Eq. (2-17) gives the null relations L°- L =io,A(K°-K ) and 


L=i o,A K hence L* = i o,A Ko. Eq. (2-23) becomes one for Ke only: 














eee eee ot En 
iAo = ido | 
1] M io, A | 4 ge 


Eq. (2-13) has to be used to find the difference electric current: 


—> + —_ 


K’ - K’ =io,[N + iMo,A|K - 2i0,H* 


Therefore, 


Zz 


(2-24) 


(2-25) 


at + i0,,N + iMo,A|\K = io,Hiy (2-26) 


=o INC 


Since the last term in Eq. (2-26) is K , 


K?* = -{1 + io,| N + iMo,A| LK (2-27) 


L* and L’ can be obtained algebraically by multiplying io, A to K° and K respectively. 


On the other hand, by Eq. (2-13), 


1 = = i, {A + |M - iNo,A|\K (2-28) 


Note that the Z = 0 case includes the special situation Z* = Z~ = Z = 4 = O when 
both sides of S are perfectly conducting. Under this circumstance L = L~ = O and the 


operator N is never involved. 


When Z is rank 1, Z =O but det Z=0. The right hand side of Eq. (2-17) provides one 


linear relation between the components of L - i o,AK which can be used to reduce the four 


components of K and L as the unknown quantities in Eq. (2-23) to three so that the 
remaining three components can be solved. The left hand side of Eq. (2-17) assures that the 
same linear relationship between the components of L - io, AK exists between 
corresponding components of the difference currents. Eq. (2-13) again has to be invoked to 


compute three other linearly independent combinations of the components of the difference 


currents from the sum currents. 
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E. CONSIDERATIONS FOR A CLOSED SURFACE 
When S is a closed surface, the choice of Z can be arbitrary for an exterior problem 


such as scattering while the choice of Z* is arbitrary for an interior problem. It is desirable 
to choose Z" = -Z* so that Z=0 and A = Z* = -Z. Then we have L = id, & K and only K 


has to be computed. With such a choice, for an exterior problem, 


ere SC ] az 5 74 
K** = ae = 0, GZ) ons (2-29) 


he 10 p ; eS " 
jo ee +iNo,Z* - M)K (2-30) 


and for an interior problem: 


— 


JR = ( = OmviaOnZe = io,N)K (2-31) 


— 


SC 10, be : ee 
DO“ = —2(M - 2 + iNo,Z-)K (23m 
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Il. A THEOREM OF ANISOTROPIC ABSORBERS 


A. AXIAL RADIATION FROM A SURFACE OF 90° ROTATIONAL 
SYMMETRY 


The xy-plane cross section of a surface § having a 90 “rotational symmetry around 


the z-axis is shown in Figure 3-1. Because of this symmetry, $ can be decomposed into four 


non-overlapping congruent pieces S,, S,, S;, S, so that a 90 °rotation will bring S; into S,, ;. 





(These subscripts are considered as equal under modulo 4.) Therefore each piece S. of S can 


15 


le ole 
be parametrized in the same orthonormal coordinates (u, v), with #@=—, V= a the 
u v 


orthonormal basis vectors on S., as follows: 


ry = (x4, v), y(u, Vv), (4, Vv) ) 
To = (aye) ee), 2 te 
(3-1) 


r, = (-x,(u,v), ~y,(u,v), Z,(u,V) ) 


r = (y (u,v), ses Ue). Z (u,v) ) 


where 7. € S,. As an example of a possible choice, u = constant and v = constant can be the 


lines of curvature of S). 


In terms of the coordinates (u, v), the sum surface current distributions on S, are: 


K(F,) = K,(u,v) 





cae s (3-2) 
Ee eV) 
Since for r » ae 
; ik|F-* | 
G(r) = Ga) 
Ko 
VG) = ikrG = -V G@ Be 


the radiation from such current distributions at a distance r » max | 7,| along the positive 


Z-axis is, from Eq. (2-6): 
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EB (2-8 oes 
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4 3-5 
_ ik ff. > [K, (u,v) a L.(u, v)] (525) 
“=I 


4tr 








4 alae 
*y), [K; (u,v) “Lupe ZW 20) *%5 *o dy dy 
1=] 


Note that this approximation is independent of the wavelength; it is applicable in regions 
closer to S than the usual Fresnel zone. 
B. CONDITION FOR VANISHING ON-AXIS BACKSCATTERING 

Consider two situations when a linearly polarized plane wave of unit strength is 
incident on S along the z-axis from the positive direction: Situation 1, identified with the 
superscript (1) has the wave polarized in the x-direction while Situation 2, identified with the 


superscript (2), has the wave polarized in the y-direction. The incident fields are respectively: 


—; (1) Za -~ik 
=~ inc (1) : (3-6) 
He y ce 

rs inc, (2) ~  -ik 

+ inc, (2) ’ (3-7) 
ee eter tc 


Note that, as seen from the positive z-axis, the incident wave in Situation 2 is that of 


Situation | rotated by 90° counterclockwise. Furthermore, Situation 2 can be obtained from 
Situation 1 through the duality transformation E'“ +H", H'°--E"°. Therefore, for 
a plane wave, a90° rotation 1s equivalent to undergoing the duality transform. 


Because of the rotational symmetry of S, the currents excited on S. under Situation 


| must appear on §,,, under Situation 2. Therefore: 


K (u,v) = - Kv) 


i+1.x y 
(2) (1) 

Kie1yv) - K; (u,v) (3-8) 
(2) (1) 

Kj. 1,2(4,) = Kj (u,v) 


Ney) = - L;, (u,v) 
LG) ec) (3-9) 
Lee v) = i (u, v) 


Assume that Z on S is invertible, the sum currents on § are determined by Eq. (2-22). The 
tangential components of the incident fields which appear on the right-hand-side of that 


equation under these two Situations are: 








ur xX 
= inc.(1) oan 
Fan a Sika 
See ae (3-10) 
’ —U° 
ae 4 
Neos 
iy 
pire ane pire 
tan Vegy _ik Oo = tan 
. = | Moe C iam | (3-11) 
Be “rx I O Naps 
tan tan 
o£ a 


where / is the 2x2 identity matrix. Therefore, 
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uM -N| le Ro Bie 
N Mi {pO nas (| ey (1) on 
iL fe jel 
uM -M\E® Fae Bee 0 -1 ele 
N M\lr@2 es SO) | wee O | | inc.) oer) 
L 6 idee 13 aes 
Since 
O -/i'|\M -N NE IQ) 7 
= (3-14) 
I O||N M N Mir O 
it follows that if 
Sil Oe 
= R (3-15) 
I O 1 «0 
ee ie 
we can multiply : to Eq. (3-12) to get: 
= (1) =(1) Ee 
Via | O Sf 0 -I\\|K O -TI} \ftan 
nv Mitr oflrml” “tr ollrmt> *lr o| fqn ~ 
18 ib Che 
Therefore the excited surface currents in these two situations are related by: 
rae 9 -1iz® AL 
= = (3-17) 
C2 I 0} lz@ Faw 





Combining this result with Eqs. (3-8) and (3-9), we have: 


(2) = (1) (1) 
EC, — a Di.) (U,V) a Ki (u,v) 


(2) = (1) 7) 
Kj yu ¥) ~— 1a Lis 1,y(UsV) cas i (u,v) 


(3-18) 
Hee (, V) = CD) — Ly (v) 


(2) — pil) = «7 Oh) 
Lis i,y(4sv) = K;.1,y (u,v) =o, en) 


so that 
4 ; Bs 
SO cep) +78) (a) = 0 (3-19) 
i=l 


and 


4 


> (Kk, @v) - Lv) = 0 (3-20) 


i=] 


Hence, along the positive z-axis, from Eq (3-5), 


. 4 
Br = 2% i A? [Ky (wv) + Luv) 
t=] 


Atr 
4 


; ae yeas a (3-21) 
+ gy [KW = Lu] fe KY Zo)" *%0 “Yo dy 


y 
t=] 


and the backscattering from S along the positive z-direction must vanish if Eq. (3-15) is 


Satisfied. 
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c: IMPEDANCE MATRICES FOR ZERO ON-AXIS BACKSCATTERING 





Z-hZ = Zor ah 
It can be verified that the matrix commutes with 
=10, Z (A Geaae 0 
if and only if: 
OSA (SAIN Zi (3-22) 
2-6, 2. Oem NZ (3-23) 


where both Z and A are 2 x 2 matrices. Under the assumption that the inverse of Z exists, 
we analyze Eqs. (3-22) and (3-23) as follows: 


Because of the identity: 


] 
detZ 


at 





Dy PAU Kale (3-24) 


and, by multiplying o,, to both sides of Eq. (3-22): 


FBI (0), NPA EC (3-25) 


Fg. (3-23) can now be transformed to ~ 


ZL, 


SSG, AcZ 50.) One rer 


= - Aa, Aa, (0,Z10,) + 6,Z 6, (3-26) 


] 
=—— (] - (Aca) a 
eo ( 20 | 
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where 


(Aa,) = Ao,Ag, = (3.2%) 


It is observed that Eq. (3-27) is greatly simplified if A 1s symmetric. Therefore, in 


this thesis, we consider only the case when A = A’ . Then Ns = ae so that: 





(Ao,y = (detA)I (3-28) 
From Eq. (3-26), 
J) = lea Y 8 
Le 3-2 
| detZ | Cn 
Eq. (3-29) can be satisfied only if | eal = +1. These two cases are: 
et 
Case I 
0 24 
Z=-Z'= » 2, #0 (3-30) 
2.6 a0 
detA = 1+ detZ =1 +z, (3-31) 
Case I 
Ze (3-32) 
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det Zo det i= al (3-33) 


On the other hand, substituting Eq. (3-24) into Eq. (3-22) yields: 


A ip ot) (aig Sy) (3-34) 
Z9(A2- Aas) = (212-29) An 6235) 
24;Ao9 + Aq, = 2%,A,, = 22,,4,, (3-36) 


For Case 1, Eqs (3-34) and (3-35) require that A,, = A,, = 0, Eq. (3-36) requires 
that A,, = A,, = 0.Therefore A = 0.FromEq (3-31), 1 + Zi, = 0. Therefore Zo = Fi 


tautinat Z = Z =Z=+ Gs 


For Case I, Eqs. (3-34) and (3-35) are trivially satisfied. Eq (3-36) becomes: 


Z4, Ao) + ZyAy, ~ 22%. 4), = 9 C7) 


Eq. (3-33) 1s explicitly: 


£11%22 ~ ae ag es - ie =). —(3-38) 


The sum of Eq. (3-37) with Eq.(3-38) is: 


(Zt Ae iz Ao) — (2 a) = det(Z alec 
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Subtracting Eq. (3-37) from Eq. (3-38): 


a Nea) (Zou Digg) az Noe = det(Z - A) =detZ =1 (3-40) 


In summary, two sufficient conditions to satisfy Eqs. (3-22) and (3-23) have been 


deduced under the assumptions that A is symmetric and Z 1s invertible. The first one is: 


Z =Z ==£40, (3-41) 


The other, with both Z* and Z symmetric, is: 


detZ == detZ> = 1 (3-42) 


It should be noted that if S$ is a closed surface, the impedance boundary condition 
closes off the interior of the surface from its exterior. Therefore for the exterior problem, 
det Z* = lis sufficient to eliminate the on-axis backscattering from a body of 90° rotational 
symmetry. Furthermore, Z* may vary with location. This is an extension of Weston’'s theory 


of isotropic absorbers [4]. 
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IV. SCATTERING OF AN ANISOTROPICALLY COATED CYLINDER 


In this chapter, the sum-difference surface current formulation developed in Chapter 
II is applied to the problem of scattering of an anisotropically coated tubular circular cylinder 
of finite length and negligible wall thickness. Due to the rotational symmetry, a Fourier series 
expansion can be utilized to reduce the variables of the problem to only the one along the 


axis of symmetry, chosen as the z-axis. The Fourier components M_, N, of the operators M 


and N are deduced in terms of the Fourier components of G(7-7,)and its partial derivatives. 
To solve the set of integrodifferential equations so obtained, the equations and the surface 
currents are both weighted and expanded over the Chebyshev polynomials. The resultant 
infinite system of linear equations are then truncated and inverted numerically. 
A. GEOMETRY AND COORDINATE SYSTEM 

Figure 4-1 shows the geometry and coordinate system of a tubular circular cylinder 
of radius a and length 2h. The cylindrical coordinate system (—, @, z) 1s scaled so that the 
surface of the cylinder S is specified by p =1 and -1 < z <1. The radial vector is thus given 


by: 


™Y 


oe 


app + hzz 


ap (cosbt + sinds) + hzé Cel, 


To facilitate representing the tangent vectors over S in matrix form, we chose i = @ and 


aA 


Vv =Z as the orthonormal basis vectors on 5S, so that 7 = f is the outward normal to S. 

To properly classify the polarization of the incident wave, the rectangular coordinate 
system (x,y,z) 1s determined as follows: When a plane wave is not incident along the z-axis, 
the coordinates are chosen so that the wave is propagating in the zx-plane incident from the 
half-plane containing the positive x-axis. Axial incidences then are considered as the limiting 
cases as the direction of incidence approaches the positive or the negative z-axis from this 


half-plane. Therefore, even for axial incidence, a linearly polarized incident wave having its 


— 


pss, 





Figure 4-1. The Geometry and Coordinate System. 


electric field intensity vector E pointing in the y-direction is considered a TE wave while one 
having its E vector in the zx-plane is considered a TM wave. In the spherical coordinate 
system (7r,0,), the direction of propagation of the incident wave k is given by 
k = -Zcos 6. - <sin®. where the incident angle 0, varies over the range 0 < O. < 7. For 


an incident plane wave of unit strength, the fields are: 


TE - polariation: 


tker 


Ine 
E = 


e 
EN A HAN Lies ole ast NG) a 


TM - polarization: 
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= Inc “Fy inc x 2 Por pee 
E™ = H™ xk = -(&cos@. - Zsin6,) e**” 


where k = kk. Note that ke 7 = [,zcos0,+1, psin8.cos@ where J, =kh, and Ll, =ka. 


On the surfaces S, the tangential components of TE-polarized plane wave are: 


E, = cose’ 
H,~ = - cosO,sind e*" (2215) 
BE = Saijal 


ré 


and the tangential components of TM-polarized plane wave are: 


E, = cos@,singe"’ 
E,™ = sin8.e""" (4-2) 


Hy” = coshe™*’” 


it 


B. SCATTERED FIELDS 
From Eq. (2-6), the scattered fields from the cylinder are given in terms of the sum 


equivalent electric current K and the sum magnetic current L: 


SEE“ =- —V x ['[* Ler) GF-¥,)db,de, 
LL, k J -o 


+i [. i K(r, )G(F-7,)d dz, (4-3) 


l lfr2an ae 
an Ny, © a 
= [ I, K(r,)°V_ G(r-r,)d,dz, 


and 


2) 


Te = <V x [ f° Rr.) Gtr F db, fm 
+ [ ine LG? G(r-r,)do dz, (4.4) 


l ] In = 
arin. L e V G(r-r.)dbdd 
k2 IL I. (r,) * V.G(r-F )db,dz, 


Their components in cylindrical coordinates are: 
4% 5c i ee 1 2m Re A = ~ = 
71, Fe O= 7 ac], [yp Orb Ho) GFF) db, de, 
] 6, ] 2 = a5 
ip 36 ier: LF) G(r-r,) db, dz, : 
tif) [* (B9O,) KF) GF-F,) do, de, (4-5) 


2n 
7 atl i K, (r,) G(r-r,) do, dz, 





aes ] 2n - -= 
K G — d d 
l is Sah Ih nti) (r r) d, G, 








4u SC Dae oO re = 
i’ ™ H eG i (6*,) Ly (F,) G(F-7,) do, dz, 
1 oO ] 20 = ws 
> hs I, I, L (7,) G(r-F,) dg, dz, 
+1 {foe b,) Ky(F,) G(F- 7) db, dz, aan 
a ™ Ky(F,) G(F-7,) db, dz, a 
i 0 zoel I, 
l l 2n 2 oe 
iL,P - of I, K(7,) G(r-r,) do, dz, 
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atlas 27 a 
i Oe ae me [Ord ) Le F,) GF -F,) dd, dz, 


rn 1 36 I [@- ,) Ly(7,) GF -F,) db, dz, 





1 f2n = aS 4-7) 
ae Ve G(r-r.)dod.d ( 
] 7 = i i o (7) (7 r) , Lo 


FPS f [KG GG-F,) dd dz 
i; come 0 ‘eae 


sc fret coe ] 2 tee =~ een: 
ee *) I aa Ir (O° ,) Ky) G-7,) db, dz, 


1°2 l 
| 0 1 21 = —- = 
Lp oo iz iL K 7) Gr i) dh, dz, 
tif’ [* (eb) LyF,) GE-F,) dd, de, (4-8) 


We 9. aye eo db d 
zap Fh (F,) G(F-F.) db, dz, 


l 


Sulla ] Uh. *7)G(F-F)d d 
ib me iz I ae >) (7 of , Zo 





Oa as [°°b) Ky F) GEF-F,) dd, dz, 
7 on [. fr K.(7,) G(r-7,) do, dz, 

-1 JO 
“if fe @ebaty (F,) GF-7,) db, de, (4-9) 
2 -~ = 
eal L,(7,) G(r -7,) dd, dz, 


ee 1 “cb = G > = d d 
LLp abaz —fif BZ). (r i) b, Lo 
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sc 1 fe) ] Moye ~ A s 
a: Carers -$,) Ky (F,) GF-F,) db, d 
We Zz (7) he mo i p (d h,) Cae, G r) , Z, 


2 2m > = 
ws x6 I f° O24) Ko G,) GF -7,) db, dz, 





l 20 x ies 4-10) 
je L db d ( 
l ° a I ih ol.) G(r-r,) dh, dz, 


“ted P ae ae [LG GG- r) dh, dz, 


Note that in the above equations, 


6 = p,cos(p-,) + o,sin(d-,) (4-11) 
@ = -6,sin(>-o,) + $,cos(b-,) (4-12) 


Because of the rotational symmetry G(r-r,) depends on @ - @, and Fourier series 


can be introduced to eliminate the variable @. Define the Fourier series expansion of a 


function f(@,z) by: 





f(o,2 = >, e™ F(z) (4-13) 
then = 
or ara co HO. 
G(F-F,)=———-= )o ee" "G1, |z-z,|. 4» 2 (4-14) 
k\r=r,| ns@ 
and 


30 


x ad -in(® - a 
Gl 12-214, 9) = f Ae 0-90) GF-F) 
-% 


Eqs. (4-5) to (4-10) become: 


SC ] 0 ] 
a 7) = ve 2 Ga + Galdz 
ne oni P z) ay az J - ll ol n=l n : 0 


~ =P L_(Z,) Gaz, 
ia p J-1 




















Fa 
+s [. Kee, Gamer Gamma, = mal K,,(z,) G,4z, 
bi ee 
ie cae a 
TL, aps faz, | Ono 
2 SC pe x l 0 ] 
Be ade @ jdz,+ —— ie G_dz 
BL. oni P z) Dy en  HonZo dt n-1 cn 1 Xo L, dp -] EA) 
+1 on(Zo) [5 6. dz, 
i ar 
-—4_ " Qe 2 
fo Ay Oz, wn. 0 "i 0 
El P.2) = T— fp) Lyle MO-VG,., -(0+ DG, dz 
iL. re 21,p _ on\*o n=) n 0 
1 Q 
2, Bp [. eG HGerlaz 
n oO 
ee Gd 
ae Oz i tn{Zo) no 
a Mou tif’ K,(@,) G,dz 
We Oz ey oz, wn. 0 +43 0 =I wa. Oo ft oO 








Sil 


(4-15) 


(4-16) 


(4-17) 


(4-18) 











1 
Bee p67) = —-—— —— G_.+G_..Jd 
LL, Sri OS) a 82 (20) Coa Gea. 
‘ in K 2) 6,42, 
Ps 
l fri re) (4-19) 
ay f, 4) [Grea G,,.,)4z, saree Ly, (Z,) G,,dz, 
= <L _,(Z,)| G,dz, 
ee dpJ- 
p SC EL oO 1 C l 
Saar). == G dz = = KK GG 
jie) =~ oF Sf, KoMGr-1~Gratldeo- 7 5 [, Kuo) G42, 
pal ] n? 
se i- Lyn(Z 4) Hien "Gua Cale (4-20) 
[ Gbal Gude 
[ip 4-1 oz 
—HX(p,2) =- <1 [! Ky(c,)-1G,_.-(n+1)G, 1) 
LL, Dy ej ‘ ; 
| axe) 1 
¥ 21, dp ie Ko(Z)(G,,-; +G,,,)dz, 
ee (4-21) 
ia foe Gd 
hb =e ie n\Zo) nAXo 


= 
2 es =f ist, CIS orgs ie LaG@o G da 


Note that in Eqs. (4-16) to (4-21),G, stands for G_(/ 





/,,@). In shifting the z-derivative 


to the current densities K, (z,) and L, (2,), the fact that edge conditions [5] require that these 


components of the scattering currents vanish as | z,| approach | from below is used 


Be 


As p ~ 1*, the operators M, and N_ can be deduced from Eqs. (4-16) to (4-21): 


2 
=“arghlal dz, 1-1 * Gus) -—G, 
[, M2 


2 

















6, - 
Myao=n f dz, [G a az, 
(4-22) 
M,,.21 = n< [i de, LG Jp =1 
1 [6 re) 
i = ~ ilyl, [ (Galen i oe alo =15 
l 
a = Ze = fi az ea C. p=l 
l, 
N a! p dz, 5 1Gale- 7“ 
a (4-23) 
N21 = 5 "dz, | G,o4-O DO rly Se: s(n “Gwe 
; 7 
N22 = 9 
Note that: 
1 oe] fe] 
ff. 21,0. CC aa aap cia 1G | seg | ea 
2 


Assuming that Z is invertible, from Eq. (2-23), the equations for the sum currents are 








therefore: a 
M, “Np K, K, Ee 
ariel see ee (4-25) 
N, M,\|L bg a 
where 
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NT a) SITS ZB atoyn 
Re (4-26) 
age 2 ARMOR 6x0, 


5 z= = INC IC : : é : 
amelithe VECtOIsh fd 3 E. Lon n ae two dimensional column matrix representations of 


the respective tangential vector fields on S over the orthonomal basis o, 2. The incident 


fields are, for TE and TM polarizations respectively: 


TE: 
inc ley : -il, zcos@; 
Ey > Gap oC esm.O)je 
, n(-i)"*cosO. eee = 
Hy, = a J_ (1, sin8,) e ead (4-27) 
, sin8. 
Hi“ = -sin@.(-i"-'J (1, sind) e 
TM: 
n(-1)"cos8@. Jee) 
i BG | ~J (1, sin6,) e cs 
l, sinO, 
Ei" = sinO.(-i)"J, (1, sin®,) e 1°"! (4-28) 
He. = (nye (Lysine) ame 


©: TRANSFORM TO SYSTEM OF LINEAR EQUATIONS 


Since the surface current components K4(9,z,) =O(1 -z) 12 and 


K(4,z,) =O0 -z" yi aaas |z,| 1”, representations of Ky, (Z,) and“ K an(Z,) in Chebyshev 
Lo 


polynomials of the first kind combined with the weighting factor conform to the proper edge 


behavior of the currents: 
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K, (Z,) = ke Pz) 
onl Tt sinv im bn “ph 





> Ki, cosqv (4-29) 
Tt Sinv g=0 


= 5 K", sin(q+1)v 


K, (2) 
we) = =O 


i 


where T (Z,)is the Chebyshev polynomials of the first kind and z) =cosv, -l<z)<1. 


Similarly, 





Ly, (20) = meres — Yb, cosqv ; 
: - (4-30) 
L, 0) = — Lyng, sin(q+1)v 


For an invertible Z, the coefficients in the above equations are to be determined from 
Eq. (4-23). To make use of the orthogonal property of the Chebyshev polynomials, the factor 
sinv sin(p + 1)vfor p>0 is multiplied to both sides of Eq. (4-23) before an integration over 
the range of v is carried out. The results are described term by term in the subsection to 
follow. This procedure creates an infinite system of linear equations to be solved numerically 


after it is truncated at an appropriate order determined by the electrical size of the cylinder. 


1. Incident Fields 


pe = Inc = 
tan,n nav . 
= ara — sinv sin(p+1)v| (4-31) 
ae Pp +] Jo TT -7 inc 
tan,n tan, 7 


a5 


TE: 


Be) Ge) ell Senco) seme (cos) | -I,(Upsinig) 


inc, 2 “e ; 
EU SL eos) i sine 
fn . a | moon sD I (4-32) 
| Losing) a 
ee = (27) ein 6; ( J,-,(1,cos 8;) + J, , (1, cos 6,)) J_(i,sin8,) 
TM: 
inc, 2n ( S 1)" fs : 
BS ee a ees ay ci 
n L, L sin 0. Ak 1 ) A 2 i) 


eee = (2) sings ( J,-(/,c0s 8;) + J, .,(2,cos 6;)) J (L,sin@.) 


(4-33) 
He, = (07 2" (eeicus3,) sey cosa) er @r sings . 


where J‘(«)is a derivative with respect to the argument. Note that as Q. approaches 0 or 7, 


only n = +1 terms are nonzero. Hence only n = +1 currents exist. For axial incident when 


B= 00): 


t 








iE: 
Inc,p  __ (au _ inc, 
1 ~ iw Joh) a Bae 
1 
Ro: (43) 
Inc,p (=r) = inc, 
Hy? = - 7 = lel 
TM: 
inc,p  _ (-ip*! 2S inc, 
Be Jay = = a 
(4-35) 





icp: = (Sn? = inc, 
eae - ] Je) 7 ey 


36 


2 The R - Matrix Term 


— [ SY sinv sin(p+1)vK,,(cos v) 
ptl/Jo tm 
2 - ae 
Sees. K;, { “dv sin(p+1)v cosqv 
(p+1)%? 2 : I 
- Dag" Ke, 
q=0 | 
where 
0 ptq odd 
Pq _ 
Aes noe Ppt+q even 
t*(p+q+1) (p-q+l) 
2 prdv . , 
—, = sinv sin(yt+l)v K, (cosv) 
- em Ki, [ "dv sinv sinp+ lv sin(q+1)v 
$i ¢=0 
eal Ks 
WE 
where 
0 pt+q odd 
Psd _ = 
A, = ee 2 p+q even 
ape tl (pot ly(ptg+3s)(p-q—)) 
Therefore: 
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(4-36) 


(4-37) 


(4-38) 


(4-39) 


Pea 
Ars 0 OmanO 




















oe K| -« ONO 
= sinv sin(p+l)v R = R 
pel L| 47-0 0 0 eA? OnE. 
(4-40) 
ee OE 
3 K? 
= > RPeq| 
q=0 in 
where R?*? is a four-by-four matrix, given in terms of R by 
Rj" = R;, Ag” 
Ree =e 
Step | uolOtiat oleae (4-41) 
Riz" = Rj; Ag 
Ri," = Ri, Ae 
Note that R;" = 0 if pt+q is odd. 
3: The M_, N, Operators 
Define the 4x4 matrix X?”" by: 
dy M, -N\\K,| = - K# 
— sinv sin(p+l)v = ee (4-42) 
iat es —— sinv sin(p+1)v M, |, [K, pad 
- ~ily I n av l n q 
= > ——— amp eon —— [Cospvcos —(G,_1+G,.)) ==. Gee 
q=0 ptl Jo - ly | 
= pee, Kg, 


q=0 
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and, 


A I : +2, cae +2 
=e 5 6Gn-t~ Gat Cem 2 2 (Gi Ca (4-43) 
2 





Pq _ 
Xn] z= 





le — sinv sin(p+1)v M, ,, [Ky Ae 
-2n nav 
q= =0 p+l on 


= pa an Ie A xpety [2c COS QV | GK, 
a ; 


—sin(p 








_ Pq 
© Xn 21 Kes 
q=0 





and, 
Xe = 2nGy'4 ik: 
pile av Sinv sin(p+l)v N rll [Ky me! 
= x dv a 
; =0 pel is aon ; ey 
ae 
= y xr Se 
qg=0 
and, 
X31 al AGH Clues | oe) 


au 


pile ay sinv sin(p+l)v N n21 [Ky nZo)d 





1 *& sinvsinpty [7 we cosey, 
q-0 Pp” 


l, 
@-l\G -@+1)G | "2 ay Got Gio 





oo 


P»4 q 
»S Xn 4l se 


t| 














q=0 
and, 
iE 
X08 = Say Ut GI GLY Om GEG 
ty me Po _ (xP *2.4 ae) 
3 =G, +\GiGas Ge ) 
pai ood es — sinv sin(+l1)v M,,, tae 
= ye ix © sinvsin(p + Ee 
a yy pl 
a Be, 
q=0 
and, 


xO = BANCE Gere Gig 


peal (4-47) 
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pe wy sinv sin(p+1)v M, ., [K, ,(z,)] 





























p+l-o 
-2il 1 av 
= y bas Fie sinvsin(p +1)v : 
q=0 pti 0 © 
dv 
= [7 sing +I = one (q+1)cos(q+1)v,G,\ K/7, 
Guar OvJo T , 
=i d, Xn. ie 
q= 
and, 
ill 4(q+1 
a ie +2,q+2 +2 +29) _ 4(Qt 1) aprlarl| . 
Gee a 5 —(Gh4+G? ad -GPa -G? 4) _ G? q (4-48) 
re iy 
i Ths — sinv sin(p+1)v Ny [K,,,(z,)] 
[ocr 
= 2 pd sinvsin@ Dy [7 —sinv,sin(q+1)v, : ae ie 
q=0 “s ot, | 
7 », X32 KY, 
g= 
and, 
Pq _ 4, Crd a ie eee ead 7 
Xx; 32 i A(p +1) 3,6 Ge =Ge wG. ] (4-49) 


Note that, for n=0, the expressions simplify to: 


4] 


il, 1, 





; ; +2. 
xP4 =~ GPt_ GPa 
pti 
P»q L P.q p+2q l, d P»q p+2.q 
X41 eer te Cy )+———(G, °-G, J 
fear | Zor 
ees, Balls 1 (Gi Gh eae Gn! 2 Gk eC ae (4-50) 
We |p eel fe 
1 
Lb @ +2,q+2 +2 +2 
xd = —(GP7+G? »q 2Gt4 2G a 
0.52 4(p +1) dl, oO O ré) re) 
Xen zs Xo i Xo 44 = Xo = 0 
Also note that X”’%=0 if p+q is odd. 
M,, oie 
From the symmetry of , we can deduce that: 
ie M, 
Kee z - X73 ee a Xn 
XO a Ps Xn 7 Xia a) (4-51) 
Xi 33 = Xi Xn 34 a X12 
Xn = xP Xi, 44 zs a 
Note that for p+gq odd, 
Ae = nat = nt = Ana X25) = ees en eee ee (4-52) 
and for p+q even, 
Xn s ees a ae 7 Xia = Xi 3 a X34 = (4-38) 


The integrodifferential Eq. (4-25) is thus transformed into the infinite system of linear 
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equations of the unknown coefficients 


RK q pimee 


> (XP4+RP9]) “l=2] 0" (4-54) 
q=0 L? Ele 


which is to be solved numerically. 
D. RADIATION IN THE FAR FIELD 
In the far field, we can write Eq. (4-3) as: 


a eee, [ ‘dz, i “do [L(F,)xF]G(F-F,) 
LL, “1 °So 


+if'de,f db, {KF -if [K,(7,)sinO sin(-,) Ce) 


+ K(F,)cosO]}G(F-7,) 


or equivalently, 
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E*(n) = il, ] on AD, 6 : - . ; 
GH 5 [of =| [K,cos Osin(h b,)-K,sin 0 L,cos(® p,)] 


can 


+ b[Kycos(-$,)-L,cos0 sin(-,)+L,sinO} } e “Tis "banded. 


i -i)"e ing Pag jyp*n 


ue 
Doe 
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nat 








iJ (1, sin 8) K, -J(Lsin®@)L? | J (1,cos8) (4-56) 
dn n\"2 dn] “p"l 


‘9 (1,sin®)sin@ K.” cy (Z,cos@) + J, (2, cos6)) 


ncotO , 








ee sin@) Ge a) (ESiml te} (J, cos 9) 
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As @-0, only the n=+1 terms are nonzero in this limit, and 


l, 


- = Jq(l,sin 8) sin 8L,, (J, (1, c088) =a J,,.2(1,cos 8) | 





06=6 =fcosh+ sind, h = -*Lsind + Pcosh. We have 6 + id = (< + S)e7?: 


Be ie 3 (-i? {2 (KZ, - Kb) +i(g, “Lg 4) 





(4-57) 
: is|(KE, +Kg.)+ (LE, Lé.)]} JM) 


A 


Similarly, as 8 — 7, only the n = +1 terms are nonzero. 8 = -6 = -fcosd - Fsind, 


aw 


@ = -£sind + f~cosd. We have 6 + ih = -(% = S)e*?: 
y 





GA 4 pa0 | | | (4-58) 
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E. INSIDE AND OUTSIDE SURFACE CURRENTS 
From Eq. (2-20) we have: 


| alien) als WEEN ee | 
|B es Oe “i0(Z-AZ Ac Aa er ig 


Therefore the following matrix equation is obtained: 


Koa Zia Zenon cr 
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(4-59) 


(4-60) 





V. COMPUTATION AND RESULTS 


The solution to the problem of the scattering of an anisotropically coated tubular 
cylinder of finite length as formulated in the previous chapter has been coded in FORTRAN 
and tested. The program listings are included in the Appendix. Computation has been carried 
out on the 32-bit Sun SPARC Station running under the Unix operating system in the 
Electrical and Computer Engineering Department and the Computer Center. The evaluation 
of the double series expansion coefficients of the Green's function and its derivatives for 
greater values of ka and kh have also been done on the 64-bit Cray Y-MP EL98 in the 
Visualization Lab so that the accuracy of the results can be accessed. | 

The program accepts the geometry of the cylinder, the surface impedance matrices, 
the incident wave frequency, direction and its polarization in terms of TE, TM or some linear 
combination of the two. It computes the sum surface currents and the far field radiation in 
the directions specified, including the strength and the phase of each field component. It also 
breaks down the sum surface currents into the outside and the inside Hess. 

In this chapter, some interesting results of computation for the scattering of a tubular 
cylinder having the length-to-diameter ratio h/a of 4 or 6 are presented. All figures are 
attached at the end of this chapter. The wave is incident from the positive z-axis and is 
polarized in the ¥ direction. 

A. COMPARISON WITH EXPERIMENTAL DATA 

For backscattering from a perfectly conducting tubular cylinder of a y-polarized plane 
wave incident from the positive z-axis, experimental] data are available [6] over a frequency 
band of well beyond three octaves, with the circumference-to-wavelensth ratio 
ka = 21 a/i varied from 0.9448 to 3.3152. These data are measured with two sets of four 
cylinders each; one set having h/a = 4, the other with h/a = 6. Both data sets use the inner 
radii of the tubular cylinders as the parameter a [7]. The cylinder length-to-wavelength ratio 


2h/X varies from 1.2030 to 4.2210 for the h/a = 4 case and from 1.8045 to 6.3315 for the h/a 
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= 6 Case. 

The experimental data are plotted against the results of theoretical computation. 
Figure 5-1 shows the backscattering from the set of cylinders having h/a = 4. Figure 5-2 
shows those data from the set of cylinders with h/a = 6. The output of theoretical 
computation follows the measured data very closely. Note that the cutoff frequency of the 
dominant circular waveguide mode, TE,,, occurs at 2h/A = 2.344 for the cylinders with h/a 
= 4 and at 2h/A = 3.516 for those with h/a = 6. 
B. NULL ON-AXIS BACKSCATTERING 

Three different ways of coating a surface impedance Z, on a tubular-cylinder of h/a 
= 6 are considered: Both on the outside surface and the inside surface; Only on the inside and 
leaving the outside surface perfectly conducting; Only on the outside and leaving the inside 
surface perfectly conducting. Here Z, has the elements z,,, = 0.5, Z,, varies from 0.1 to 5, Zz, 
= Z>, = 0. At a fixed frequency for which 2h/A = 3.194, slightly below the TE circular 
waveguide dominant mode cutoff of 3.516, the scattered fields are plotted in Figures 5-3 
through 5-5. Figure 5-3 shows the results of computation for the case Z* = Z = Z,.=Z. As 
>) 18 varied through 2, the backscattering cross section vanishes as predicted in Chapter 3. 
Figure 5-4 shows the results for the case Z* = 0 and Z = Z.. As the impedance on the inside 
surface is increased, the excited field inside the tubular cylinder is dissipated and the 
backscattered power decreases exponentially. Figure 5-5 shows the results for the case Z* = 
Z, and Z = 0. The backscattered power drops off rapidly at first as the impedance on the 
outside surface is increased. But the cross section quickly settles down to a fixed value 
presumably due to the current excited on the perfectly conducting inside surface of the 
cylinder. “ 

Results of computation for the same configurations but at a higher frequency for 
which 2h/i = 4.865, above the TE,, circular waveguide dominant mode cutoff of 3.516, are 
plotted in Figures 5-6 through 5-8. Now that the incident wave can propagate through the 


cylinder in the dominant waveguide mode, the backscattering cross sections are about an 
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order of magnitude smaller than in previous cases. Figure 5-6 shows the results of 
computation for the case Z" = Z = Z, = Z. Again the backscattering cross section vanishes 
aS Z,, is varied through 2. Figure 5-7 shows the results for the case Z* = Q and Z = Z, while 
Figure 5-8 shows the results for the case Z* = Z, and Z = 0. It appears that, above cutoff, the 
contribution to the backscattering cross section from the inside of the tubular cylinder is 
minimal: once the outside current is reduced by the increase in surface impedance, the 
backscattering cross section 1s reduce by more than 10 dB as shown in Figure 5-8. When the 
impedance coating is applied in the inside surface, the maximal reduction in the 
backscattering cross section 1s only about 1.2 dB. 
C. FREQUENCY DEPENDENCE 

The axial backscattering of the two cases when the tubular cylinder 1s coated only on 
the inside or only on the outside with z,,, = 0.5 and z,, = 2 are investigated for different 
frequencies with 2h/A varying from 0.1 to 7.5. Figure 5-9 shows the results of the case when 
only the inside surface is coated so that the backscattering is mainly due to the current 
excited on the outside surface. The reflection from the ends of the cylinder causes the 
fluctuation in backscattering cross section. Being waves in free space on the outside of the 
cylinder, the maxima and minima are evenly spaced with the minima occurring when 2h/A 
is a multiple of half integer. Figure 5-10 shows the results when only the outside surface is 
coated and the current on the inside surface dominates the contribution. The distinct feature 
in this case 1s that the backscattering cross section does not fluctuate with varying frequency 
below the waveguide mode cutoff. The incident wave is able to penetrate deeper into the 
cylinder with increasing frequency, resulting in a constantly rising strength of the 
backscattered field. Once beyond the circular waveguide mode cutoff, the wave can pass 
through the cylinder in the TE,, mode and the backscattering diminishes. The oscillation in 
the cross section at these higher frequencies represents the interference of reflected waves 
at the ends of the tube and the separation between maxima and minima should be determined 


by the guide wavelength at the particular frequency. These two situations should be 
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compared to Figure 5-11 which shows the results when both sides of the cylinder are 
perfectly conducting and the current can flow freely. The distinct notch in the cross section 
near the TE,, mode cutoff at 2h/A = 3.516 and the subsequent faster variation in the cross 
section shows the combination of the two distinct features of Figures 5-9 and 5-10. 
D. COMPUTATION ACCURACY 

The main difficulty encountered in the computation is the evaluation of G/""(/,, 1) 
and its /,-derivative by double power series sum when J, becomes large. Computations for 
Figure 5-11, 5-13 and 5-15 use the G’’? values evaluated with the Cray computer which has 
a 128-bit double precision number. Computations for Figures 5-12, 5-14 and 5-16 use the 
G** values evaluated with the Sun SPARC Station which has a 64-bit double precision 
number. For 2h/A greater than about 6.2, the SPARC Station fails to provide accurate 
results. 

Figures 5-13 and 5-15 are axial backscattering from a cylinder of h/a = 6 coated with 
the impedances having the elements z,, =Z,. = 0.5, Z° =z, =04, 2. =%, =Zp =Z)) 
= 0.3. Figure 5-13 shows the co-polarized backscattered field while Figure 5-15 shows the 


cross-polarized backscattering. 
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Figure 5-6. 
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Figure 5-7. 
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Y-component of axial backscattering (perfectly conducting cylinder , h/a = 6, Cray G) 
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Y-component of axial backscattering (4/a = 6, Cray G) 
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VI. CONCLUSIONS 


In this thesis, the sum-difference surface current formulation is introduced for solving 
electromagnetic boundary value problems when impedances are specified on both sides of 
a surface. For an impedance coated body, the body can be treated as being a surface 
separating the space into two regions of identical medium. For an exterior problem, the 
impedance normalized to the medium on the inside surface, Z , can be chosen arbitrarily; 
and for an interior problem, that on the outside surface, Z~, can be arbitrary. The choice 
when Z° = -Z~ 1s Of particular interest because the integrodifferential equation has only the 
sum of the equivalent electric surface currents on the outside and the inside surfaces as its 
unknown to be solved. 

This formulation preserves the duality nature of Maxwell’s equations and carries it 
over into the algebraic form of the integrodifferential operators in the equations for the sum 
currents. Since a 90° rotation is equivalent to undergoing a duality transform for an incident 
plane wave, this particular symmetry in the algebraic form of the operators leads to the 


sufficient conditions that if Z° = Z =+0,, or if Z° and Z™ are symmetric and 


2 
det Z~ = det Z° = 1, the on-axis backscattering of an anisotropic impedance coated scatterer 
having a 90° rotational symmetry will be eliminated. Note that in the symmetric case, Z~ 
and Z~ may vary with location. This is an extension of Weston's result [4] for which the 
surface impedance is isotropic. 

A FORTRAN program has been written which accepts the geometry of the cylinder, 
the surface impedance matrices, the incident wave frequency, direction and its polarization 
in terms of TE, TM or some linear combination of the two. It computes the sum surface 
currents and the far field radiation in the directions specified, including the strength and the 
phase of each field component. It also breaks down the sum surface currents into the outside 


and the inside currents. 


The results of computation using this program agree with measured data of 


67 


backscattering from conducting tubular cylinders over a frequency band of more than three 
octaves. For a cylinder coated with surface impedance matrices satisfying the criteria for null 
on-axis backscattering, the numerical computation also validated the theoretical assertion. 

Difficulties have been encountered about the computational accuracy in the 
evaluation of the double series Chebyshev expansion coefficients of the Green's function 
Ee ;:4,) and its l,-derivative by double power series sum when the length of the cylinder 
is large compared to the wavelength: a compiler with a 64-bit double precision number can 
only handle a cylinder having a length up to about 6.2 wavelengths. Further work to explore 


the feasibility of asymptotic evaluation of these coefficients 1s recommended. 
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APPENDIX PROGRAM LISTING 


A. INCLUDE FILES 


Sat a a a a Pt Pag Pe Pag Pt Pag Pt ag Mag Pag Pag Pag Pag Pt Peat Pag Pag Pent Binge Png Pings Pag Pag Pat Paa Punt Bing Bing Binge Pag Prag Bngp Brag Brag Bing Braap Braap Binge Bigg Bing Bina Bragg Bing Bing Bing Bing Ding Bing Bing Bing Bn Bing Bing ng 


REALTP . INC 

es TYPE STATEMENTS FOR REAL AND INTEGERS AND DEFINITIONS OF CONSTANTS. 
IMPLICIT DOUBLE PRECISION (A, 5B, eD-H, 0-2) 
IMPLICIT INTEGER*4 (1I-N) 
PARAMETER (PI=3.14159265358979323846264338327D0,PI2=PI+PI, 


+ PISO=F iar) 
PARAMETER (ONE=1.D0,TWO=2.D0,THR=3.D0,HXD=16.D0, ZERO=0.D0, 
+ DEGPI=180 .D0, EPS8=2 .220446049250313D-16) 


PARAMETER (ONEN=-ONE, HALF=ONE/TWO, THIR=ONE/THR, QUAR=HALF* HALF ) 


Pet gt gt Mgt Mgt Mgt Mngt Png Pag Pag Peat Pngp Mngg Png Pigg Pngp Mngt Png Png Pngg Pingp Pigg Png Png Pig Pingp Binge Ping Page Binge Png Page Png Png 


CMPXPT ENC 

C IMPLICIT TYPE STATEMENT AND CONSTANTS FOR COMPLEX NUMBERS. 
IMPLICIT DOUBLE COMPLEX (C) 
PARAMETER (CZERO=(0.d0,0.d0),CONE=(1.d0,0.d0)) 
PARAMETER (CONEN=(-1.d0,0.d0),CI1=(0.d0,1.d0) ,CI2=(0.d0,-1.d0)) 


at agp Pag Pag Pag Page Pag Ping Ping Page Pngp Ping Page Pingp Ping Png Page Ping Page Png Ping Prat Png Pint Binge Pant Prat Ping Page Ping Bing Ping Bing Ping Bing Bing Bina Bing Bingp Binge Bing Bing Pinay Bing Binge Prnat Prat Pant Prt Pat Pant Prt Big Bing Prat Bing Pint Bing Prat Punt 


Pet Pet Pt Pt Cet Pt Pet Pt Pag Pag Page Pag Page Page Paget Pag Pant Pag Pant Pag Png Pant Pang Page Pag Pagf 


B. INPUT DATA FILES 


gt gt Png Png Mngt Pine Prag Ping Png Bing Prat Page Png Pn Bing Pra Png Png Ping Pragf Png Png Prat nap Ping Pag Ping ng Png nge Pnge Png Png Prag Png Ping Page Ping? Ringe Binge Page Page Png Ping Fringe Page Pine Frage Pngp Pngp Pint Pint Pint Png Png Ping Binge Binge Page Pot Pre 


CLYGEOM. PRM 


30 Maximum kh (integer) 
5 Maximum ka (integer) 
8 NREGNS (integer) 


er LE LPT PRM 


1024 floating point zero IOBIT (1024 bits) 
64 floating point precision IFPBIT (64 bits) 
INPUTDAT . PRM 
.1D-1 DKAO, minimum ka (REAL) in the computation 
.875D-2 DINCA, increment of ka (REAL) 
40 NSTR, Start point (INTEGER) in the computation 
479 NEND, end point (INTEGER) in the computation 
6.D0 RHA, ratio of h and a- (REAL) 
A IE, if incident wave is TE-polarized it is 1, otherwise 0 
0 IM, if incident wave is TM-polarized it is 1, otherwise 0 
G2D0 THETAI, incident angle (REAL) is limited to 0 to 90 degree 
6 NTHTAO, total number (INTEGER) of angle of THATA 
3 NPHI, total number (INTEGER) of angle of PHIO to be computed 
0.DO THETAO, initial theta angle (REAL) 
1. DO DELTHO, increment of theta angle (REAL) 
0. DO PHIA, initial theta angle (REAL) 
cZ5D0 DELPHI, increment of phi angle (REAL) 
0 IK, set 1 to compute scattering currents on outer and inner surface 
IMPEDNCE. PRM 
0 IZ, if perfect conducting, IZ=0; otherwise, IZ=1. 
(.5D0,70200) Impedance Z+ (phi phi) 
(,4D070 2D0) Impedance Z- (phi phi) 
(-3D0, 0. 2c) Impedance Z+ (phi 2) 
(23007707 D0) Impedance Z- (phi 2) 
(2S b07 On DO) Impedance Z+ (z phi) 
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sO), 0.00) Impedance Z- (z phi) 
@a4D0;, 0. D0) Impedance Z+ (2 2) 
C25007.0 ..D0) Impedance Z- (z 2) 


ee a a a ag gt ee ag ag ag gt gt Rag at Rigs Ringe Mag Pp Rings Tagf Biggs Piaf Bigs Bag Tings Rings Pag Mags Ming Tings Bing Tings Miings Tia Meas Bing ing p ays Tangy Ming Ming Pigs Ming Tag Ming Tang Tiong Tingf Tras Mingp Mngt Ming Tiong ings Tima 


Cc. SETUP PROGRAM AND CREATED FILES 


PROGRAM SETUP 
CERREKKEKKEKKKEKEKKEKEKKEKEEEKKEKEREKRKRKKEKEKKKEEERKEKKE KEKE KEKKEKEKEEEKKEKEKKE K 
C NOTES:The format statement 1001 need to be revised if other than 
double precision real numbers are used for DKHMAX and DKAMAX. 
Statement 1001 needs to be revised if KHMAX or KAMAX exceeds 
SmeutcgEes. 
The format statement 1002 need to be revised if other than 
double precision real numbers are used for FZERO AND PRECSN. 
Statement 1002 needs to be revised if IOBIT exceeds 6 digits 
or IFPBIT exceeds 4 digits. 


KK KKK KEKE KKK KEKE KEKE KKK KKK KKK KKK KKK KEK KEKE KK KEKE KKK KEKE KEKE KKK KEKE KKK KEKE KKK KKK KK KKK 


GCGire O @Fa°a a a 


INCLUDE ‘REALTP.INC' 
INGEUDES -EMPATE. INC? 
OPEN (20, FILE='‘CYLGEOM. PRM‘ , IOSTAT=I0OS , STATUS='OLD’* ) 
TE ChOS = NEoe0)) THEN 
WRITE(”=, 3000) 
9000 FORMAT ('Cannot find the file CYLGEOM.PRM containing the ',/, 
+ ‘maximum values for ka and kh, and the parameter NREGNS. ') 
STOP 
END IF 
READ (20,*) KHMAX 
READ (20,*) KAMAX 
READ (20,%*) NREGNS 
GEOSE, (20) 
OPEN (20,FILE=‘CYLFLPT. PRM' , IOSTAT=IOS, STATUS='‘OLD' ) 
IF(IOS .NE. 0) THEN 
WRITE(*, 9001) 
9001 FORMAT ('Cannot find the file CYLFLPT.PRM containing floating',/, 
+ =" point zero birt LOBIT and precrsron IFPBIT.') 
STOP 
END IF 
READ(20,*) DTOBIT 
READ(20,*) LEPBEIT 
CLOSE (20) 
OPEN (21, FILE='LIMITS.INC',STATUS='UNKNOWN' ) 
WRITE (21,1001) KHMAX, KAMAX 
WRITE (21,1002) IOBIT, ISPBIT, IFPSiT 


ELOSE (21) 
1001 FORMAT (6X, ‘PARAMETER (DKHMAX= ‘,I3,'.D0, DKAMAX= ',I3,'.D0)') 
1002 FORMAT (6X, 'PARAMETER (FZERO=',16,'.D0O0, PRECSN=',14,°.D0.,', 
=e heaPBrm— , 14 %)et-) 
& 
C Part 2 
DKHMAX=KHMAX 
DKAMAX=KAMAX 
ONEDEL=ONE-EPS8 
KQDIM=INT ( (DKHMAX/ PI) *NREGNS+ONEDEL) 
KNDIM=INT ( (DKAMAX/ TWO) *NREGNS+ONEDEL) 
KQDIM=MAX (KQDIM, 1) 
KNDIM=MAX (KNDIM, 1) 
OPEN (21,FILE='MAINDM.INC', STATUS='UNKNOWN' ) 
WRITE (21,1003) NREGNS, KNDIM, KQDIM 
WRITE (21,1004) 
WRITE (21,1005) 
WRIE (21,1006) 
CLOSE (21) 
1003 FORMAT (6X,'PARAMETER (NREGNS=',12,', KNDIM=',I3, 
+ Poe CO) BG ee Oe ee) 
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1004 FORMAT (6X, 'PARAMETER (KNDIM1=KNDIM+1, KQDIM1=KQDIM+1) ') 
1005 FORMAT (6X, 'PARAMETER (KXCRT=2*KQDIM1, KCRNT=4*KQDIM1) ') 
1006 FORMAT (6X, 'PARAMETER (MAXNG=KNDIM1, MAXPEG=KQDIM/2+1,', 


+ ' MAXPOG=KQDIM1/2)') 
Cc 
C Part 3 
FZERO=IOBIT 
REFC=FZERO* LOG (TWO) ~LOG (PI2) -ONE 
REFH=HALF*® (REFC-LOG (DKHMAX) ) 
REFA=HALF* (REFC~-LOG (DKAMAX) ) 
C 


KO2=KQDIM+2 
DMXM=K02 
DO 100 WHILE ( (DMXM+HALF) * (LOG (DMXM/DKHMAX)+ONEN) .LT. REFH) 
DMXM=DMXM+ ONE 
100 CONTINUE 
MXMREG=INT (DMXM) 


c 
MXMSNG= INT (TWO* DKHMAX+ONEDEL) 
MXMSNG=MAX (IFPBIT, MXMSNG , KQ2) +KNDIM+1 
OPEN (21,FILE='GPQNDM. INC' , STATUS='UNKNOWN' ) 
WRITE (21,1009) MXMREG, MXMSNG 
CLOSE (21) 
1009 FORMAT (6X,'PARAMETER (MXMREG=',14,', MXMSNG=',14,')') 
STOP 
END 
GPONDM. INC 
PARAMETER (DKHMAX= 30.D0, DKAMAX= 5.D0) 
PARAMETER (FZERO= 1024.D0, PRECSN= 64.D0, IFPBIT= 64) 
BIMI TS. INC 
PARAMETER (DKHMAX= 30.D0, DKAMAX= 5.D0) 
PARAMETER (FZERO= 1024.D0, PRECSN= 64.D0, IFPBIT= 64) 
MATNDM. INC 


PARAMETER (NREGNS= 8, KNDIM= 20, KQDIM= 77) 
PARAMETER (KNDIM1=KNDIM+1, KQDIM1=KQDIM+1) 
PARAMETER (KXCRT=2*KQDIM1, KCRNT=4*KQDIM1) 
PARAMETER (MAXNG=KNDIM1, MAXPEG=KQDIM/2+1, MAXPOG=KQDIM1/2) 
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D. PROGRAM MAIN 


PROGRAM MAIN 
kkk kkkKkKkkKkKkkkKkkKeKkkKeKkKekKekkKeKKkKKKkKKkKkKkKkKkKkKkKkKkKhKhKR KK KKK K KKK KK KKK KKK KKKKKRKKKKKKKKKKKKKKKKKEKK 

INCLUDE 'REALTP.INC' 

INCLUDE *CMPATP.INC' 

COMMON /GCONST/ DKH, DKA, HH,HA,HSQ, DASQ, HSQN, DASOQN, HHSQ, HDASQ, 
+ RAH, RAHSQ, DHA, DAH 

COMMON /INPUT1/ DKAO,DINCA,NSTR, NEND, RHA 

COMMON 7INPUTZ/ IE, IM, THETAL, THESINI, THECOST, RIHEL 

COMMON /INPUT4/ 1I2Z,IK,IS,NYSM 


CALL CHKINPUT 

OPEN (21,FILE='rhzz41.dz', STATUS= 'UNKNOWN' ) 
DO 8001 IS=NSTR,NEND 

DS=IS 

DKA=DKA0+DINCA*DS 

DKH=DKA* RHA 

CALL MAXODM 

CALL RAPQMT 

CALL GNPQFN 

CALE, SCRSEFL (TR, CESTH,CESPH) 


7] 


CALL RCSPAREA (CESTH, CESPH) 
8001 CONTINUE 

GECSE (21) 

STOP 

END 


Ee. SUBROUTINE CHKINPUT 


SUBROUTINE CHKINPUT 
CER RRKKKKERKREKKEKEKEKKEKKEKEKKEKEKRKEKKEKRERKEKKKEKRKEKKEKKEKKKEKKKEKKKKKKKKKKK KKK RK KR KK KKK KK 
INCLUDE ‘REALTP.INC‘ 
INCLUDE 'LIMITS.INC' 
COMMON /INPUT1/ DKAO,DINCA,NSTR,NEND, RHA 
COMMON /INPUT2/ IE,IM, THETAIL, THESINI, THECOSI, RTHEI 
COMMON /INPUT3/ RTHE,RDELT, RPHI, RDELP,NTHTAO, NPHI, THESIN, THECOS, 
+ RHPI 
COMMON /INPUT4/ IZ,IK,IS,NYSM 
CERRKRKRKRKRKKEEKKERERKEKEKKKEKEKKEKKEKKKKKEKKEKKKEKKEKKKKEKKEKKKEKRKEKEKKKKKKKEKK 
OPEN (20, FILE=‘INPUTDAT. PRM' , IOSTAT=IOS , STATUS='OLD') 
IF(IOS .NE. 0) THEN 
WRITES(*;*) “Fail €o epen input Eile INPUTDAT.FrRE 
STOP 
END IF 
CER KEKKEKKEKEKEKR EK KKK KEK KKK KK KK KEKE KKK KEKE K KEKE KEKE KE KKEKKKKKEKEKKKEKKKKKEKKKKEKE 
C Input kh and ka. These values are passed to other 
C parts of this program through the common block /INPUT1/. 
READ(20,*) DKAO 
READ (20,*) DINCA 
READ (20,*) NSTR 
READ (20,*) NEND 
READ (20,*) RHA 
CER RRR K KER KEKE KE KHER KKK KK KKK KEKE KKK KKKEKKEKKKEKKKEKEKEKKEKKEKEKEKKEKKEKKKK KEK 
C Check against maximum kh and ka values. 
NBTW=NEND-NSTR 
DKHO=DKA0* RHA 
DKA1=DINCA* NBTW+DKAO 
DRHI—-DKAL* ROA 
IF (DKH1 .GT. DKHMAX) THEN 
IF (DKA1 .GT. DKAMAX) THEN 
WRITE(*,*) ‘Both kh and ka values exceed the maximum allowed. ' 


ELSE 

WRITE(*,*) ‘The input kh value exceeds the maximum value allowed. ' 
END IF 

WRITE(*,*) 'The execution is terminated. ' 

CLOSE (20) 


STOP 

ELSE IF (DKAl1 .GT. DKAMAX) THEN 
WRITE(*,*) ‘The input ka value exceeds the maximum value allowed. ' 
WRITE(*,*) 'The execution is terminated.’ 


CLOSE (20) 
STOP 
END IF 


Te ((DINCH .LT. ZERO) .OR. (DINCA “3ET. ZERO} een 
WRITE(*,*) ‘The increment DINCH or DINCA is less than zero.' 
WRITE(*,*) 'The execution is terminated.' 
CLOSE (20) 
STOP 
END IF 
CER EREE EERE REKREEKREKEKEKEKEKEKKEKEEEREEREE REREEERKKRKEKKEKEKK ERE EEE EE 
Input the incident angle (THETAI) and polarization (TE or TM) of the 
incident wave. The incident angle is limited to 0 to 90 
degrees. SIN(THETAI) and COS(THETAI) are computed also. These 
Parameters are passed to other parts of this program through the 
common block /INPUT2/. 
READ(20,*) IE 
READ(20,*) IM 
READ( 20.) THETAL 


2 eS Me me! 


2. 


C Check the input value of incident wave 
EF (( LE 2NE. 1) «AND. (IE. NE. 0) Phen 
WRITE (*,*) ‘Improperly specified the polarization of incident ' 
+ ‘wave. program is stoped. ' 
STOP 
ELSE IF ((IM .NE. 1) .AND. (IM .NE. 0)) THEN 
WRITE (*,*) ‘Improperly specified the polarization of incident ' 
+ ‘wave. program is stopped. ' 
CLOSE (20) 
STOP 
END IE 


? 


EP ((THETAT .LT. ZERO) .OR. (THETAIL .GT. 90.) ) THEN 
WRITE (*,*) ‘Improperly specified incident angle. ', 
+ ‘program is stopped.' 


CLOSE (20) 
STOP 
END IF 


C Calculate SIN(THETAI) and COS (THETAI) 

RTHEI=THETAI * PI/DEGPI 

THESINI=SIN (RTHEL) 

THECOSI=COS (RTHEI) 
CEKEEKEKKEKKKKKEKKEKKEKEKKKKEKKKEKKKKEKKEEKKEKKKEKKEKKEEKEKKEKKEKKKEKEKKEKEKKEKKEKKEKKEKKKKKKK 
Input the angles theta and phi at which the scattered fields are 
to be computed. They are specified in terms of the initial theta 
(THETAO) and phi (PHIO) angles, their respective increments DELTHO 
and DELPHI, and the total numbers of angles NTHTAO and NPHI to be 
computed. Thus NTHTAO and NPHI must be integers greater than 1. If 
either NTHTAO=0 or NPHI=0, no bistatically scattered fields will be 
computed. Note that the scattered electric field components are 
computed for all the phi-angles at a fixed theta, before the 
theta-angle is varied. All angles are specified in degrees. Theta is 
limited to 0 to 180 while phi is limited to 0 to 360 degrees. 
SIN(THETAO) and COS(THETAO) are computed also. These parameters are 
passed to other parts of this program through the common block /INPUT3/. 

READ(20,*) NTHTAO 
READ(20,*) NPHI 
PPeeNTHTAD LT. 0} .OR. (NPHI LT. 0)) THEN 
WRITE(*,*) ‘Improperly specified number of output angles. ', 
+ ‘Program is stopped.' 
STOP 
ELSE IF ((NTHTAO .EQ. 0) .OR. (NPHI .EQ. 0O)) THEN 
CLOSE (20) 
WRITE(*,*) ‘Desired bistatic scattered field direction has not ', 
+ ‘been (properly) specified, they will not be computed.' 
NTHTAO=0 
NPHI=0 
THETAO=ZERO 
DELTHO=ZERO 
PHIO=ZERO 
DELPHI=ZERO 
ERS E 
READ (20,*) THETAO 
READ (20,*) DELTHO 
READ(20,*) PHIA 
READ (20,*) DELPHI 
END IF 
READ(Z0*) 2K 
CEEKKREKKKKKKEKEKKEKKKKKEKEKK KEKE KKKKKKKEKKKEKKEKKKEKEEKKEKKEKEKKKEEKEKEKEKKKKK 
C Check input values. 
IF ((DKHO .LE. ZERO) .OR. (DKH1 .LE. ZERO) ) THEN 

WRITE(*,*) ‘Invalid kh, program is stopped.' 

STOP 

END IF 

IF ((DKAO .LE. ZERO) .OR. (DKA1 .LE. ZERO) ) THEN 
WRITE(*,*) ‘Invalid ka, program is stopped. ' 
STOP 

END IF 


Mo OG @ 1) OO 63 01) 


US 


IF ((DKAO .GT 


. DKHO) .OR. (DKA1 .GT. DKH1)) THEN 


WRITE(*,*) ‘ka/kh > 1, program is stopped.' 


STOP 
END IF 


C Output angle checking not required: 
IF (NTHTAO .EQ. 0) GO TO 200 
C Checking output angles: 
TE SUCTHETAO .LT. ZERO) “JOR. (LHETAGO  .CGT. DEGPI).) THEN 
WRITE(*,*) ‘The first output theta~angle lies outside the 0 to ', 
+ ‘180 degrees range. Program is stopped.' 


Se 
END IF 


THTAIF=THETAO+ (NTHTAO-1) * DELTHO 
EE ((THTAIF .LT. ZERO) @30R )(THTALE SGT. DEGPI)) THEN 


WRITE(*,*) ‘Some 
+ ‘outside the 0 
STOP 
END Le 
PHIMX=TWO* DEGPI 
IF ((PHIO .LT 
WRITE(*,*}) “he 


of the specified output theta-angles lie', 
to 180 degrees range. Program is stopped. ' 


. ZERO) .OR. (PHIO .GT. PHIMX)) THEN 
first output phi-angle lies outside the 0 to ', 


+ '360 degrees range. Program is stopped.' 


SLOP 
END IF 


PHIF=PHIO+ (NPHI-1)- DELPHE 
IF (({PHIF JLT. ZERO) @jeCh wea EHEr CT. TWO :DECPI) a THEN 


WRITE(*,*) 'Some 


of the output phi-angles lie outside', 


+ ‘the 0 to 360 degrees range. Program is stopped. ' 


STOP 
END IF 

200 CONTINUE 
DPI=PI/DEGPI 
RTHE=THETAO* DPI 
RPHI=PHIO* OPI 
RDELT=DELTHO* DPI 
RDELP=DELPHI*DPI 
RHPI=90.*DPI 
THES IN=SIN (RTHE) 
THECOS=COS (RTHE) 
RETURN 
END 


F. SUBROUTINE MAXODM 


SUBROUTINE MAXODM 


‘e KKK KKK KK KK KEKE KK 


INCLUDE 'REALTP. 
INCLUDE 'MAINDM. 


KEKE KEKE KEKE KREKKK KKK KKK KKK KEK KK KERR KEKE KKK KEKEHK 


ENG * 
INC' 


COMMON /CRNTDM/ NMAX,MXNG, IQMAX, IQMAX1, IQMAX2, IXCRNT, ICRNT,MXQEG, 


re MXQOG 
COMMON /GCONST/ 


DKH, DKA, HH, HA, HSQ, DASQ, HSQN, DASOQN, HHSQ, HDASQ, 


* RAH, RAHSQ, DHA, DAH 


COMMON /RTHETA/ 
COMMON /INPUT1/ 
COMMON /INPUT2/ 
COMMON /INPUT3/ 
+ 


ONEDEL=ONE-EPS8 


DE TCOST ,DLZSINI DOLL 

DKAO, DINCA, NSTR, NEND, RHA 

TE, IM, THETAL, THESINI, THECOST, RIBS 

RTHE, RDELT, REHI, RDELP, NTHTAO, NPHT, THEStN, Tiieeeoe 
RHPI 


IQMAX=INT ( (DKH/ PI) *NREGNS+ONEDEL) 


IQMAX=MAX (IQMAX, 
IQMAX1=IQMAX+1 
IQMAX2=IQMAX1+2 
IXCRNT=2 * TIOMAX1 
ICRNT=4* IQMAX1 
MXQEG=IQMAX/2+1 
MXQOG=IQMAX1/2 


tj 
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NMAX=INT ( (DKA/TWO) *NREGNS+ONEDEL) 


NMAX=MAX (NMAX, 1) 
MXNG=NMAX+1 


C Evaluate SIN and COS functions to pass along / 


DL1COSI=DKH* THECOSI 
DL2ZSINI=DKA* THESINI 
DLL=DKH*DL2SINI 


C Evaluate geometrical values to pass along /GCONST/. 


HH=HALF* DKH 
HA=HALF* DKA 
HSQ=DKH* DKH 
DASQ=DKA* DKA 
HSQN=-HSQ 
DASQN=~-DASQ 
HHSQ=QUAR*HSQ 
HDASQ=QUAR* DASQ 
RAH=DKA/DKH 
RAHSQ=RAH* RAH 
DHA=DKH*HA 
DAH=HH* HA 


RETURN 
END 


G. SUBROUTINE RAPQMT 


SUBROUTINE RAPQMT 


CER RK KK KEKE KKK KKK KKK KKK KK KKK KK KR KK KK KKK KKK KEK KKK KKK KK KKK KK KK KKK KK KKK K 


INCLUDE ‘REALTP.INC' 
ENGHUDE | CMPATP. INC’ 
INCLUDE 'MAINDM.INC' 
DIMENSION CO(2,2), CIR(2,2) 


DEMENS@eN GZSUM(2,2),CZDIF(2,2),CR1(4,4),CR2 (4,4) ,CR3 (474) 
DIMENSION CRPQ1 (KCRNT, KCRNT) , CRPQ2 (KCRNT,KCRNT) , 

+ CRPQ31 (KCRNT, KCRNT) , CRPQ32 (KCRNT, KCRNT) 

COMMON /INPUT2/ IE,IM, THETAI, THESINI, THECOSI, RTHEI 


COMMON /INPUT4/ IZ,IK,IS,NYSM 
COMMON =/ XPOTMP1/ CRPOL,CRPO2 
COMMON /XPQTMP2/ CRPQ31,CRPQ32 


COMMON /CRNTDM/ NMAX,MXNG, IQMAX, IQMAX1, IOQMAX2, IXCRNT, ICRNT, MXQEG, 


+ MXQ0G 
‘o 
OPEN (20, FILE='IMPEDNCE.PRM' , IOSTAT=I0S, STATUS='OLD' ) 
TF Cros -NE. 0) THEN 
WRoLeE Ge oO00) 
9000 FORMAT ('Cannot find the file IMPEDNCE1.PRM containing the ',/, 
+ ‘impedances of the inner and outer surfaces.') 
STOP 
END IF 
cc 


READUZO. * ) FZ 

ie 2CGuZ oN .. 0) 
WEL *, =) 
+ ‘surface. program is stopped.' 
CLOSE (20) 
STOP 

END IF 

TF (rz. «EO. 
CLOSE (20) 
RETURN 

END IF 


-AND. (IZ .NE. 


0) THEN 


C Set up the matrix when Z and Delta are diagonal, 


C Initialize the matrix CRPQ1 
BOr2 00 i=l KCRNT 


i) ) THEN 


1 


‘Improperly specified the impedance of cylinder ', 


or not diagonal but mec o- 


100 
200 


300 
400 


BOe100 J=1,KCRNT 
CRPQ1 (I,J) =CZERO 
CRPQ2 (I,J) =CZERO 
GCRPOST(1L.J )-CZERO 
CRPOQ32 (I,J) =CZERO 
CONTINUE 
CONTINUE 
DO 400 I=1,2 
DOr300 -J=15,2 
READS 20, je COd Ly) 
READS ZO, > | CERT.) 
GZs0Mtt 3) =(CO(l a )+CIR (1 7) eae 
CZDTR(L,. 7) =(COA1 J) -CIR(1,3)) SHALE 
CONTINUE 
CONTINUE 
CHOSE (20) 
CZUT=CZSUM(1,1) 
CZ12=CZSUM (1,2) 
GZz2l=CZ2SUM(2Z, 1) 
CZ22=CZSUM (2,2) 
GRBOET=CZ11*CZ22-CzZ12°C221 
CDLI-CZDIF (1,1) 
CbiZ2=CZDIF (1,2) 
Giza =CZDIF (2,1) 
GBZ2=C2ZDIF (2,2) 
ir (CRDET .50, CZERO) THEN 
WRITE, 4", 9001) 
SLor 
END IF 


C Check the diagonalization of the impedance matrixes 


+ 


+ 


+ 


Jt, 


NSYM=1 
TP (e212. .NE. CZERO) THEN 
NSYM=0 
ELSE IF (CZ21 .NE. CZERO) THEN 
NSYM=0 
BESE ff (€Dl2 .NE. CZERO) THEN 
NSYM=0 
BESE Ir (Cb21 . NE. CA2ERO)}s fHeN 
NSYM=0 
END IF 
CRDTZ=CONE/CRDET 
CRI (1,1) =CZ11= (CD11 *ep1 1 *Ca22-Cbl1*CD12* C221) -Cb ih Ci easy 
+CDi2* CD21 4 C711) SeCRDIZ 
CRI G2) =CZ122( CD11 * Cpl aioe 2 Cl 2*CD12*CZ21=-Cbi epZZ -eal2 
+CDI2* CD22 76201) *CRDTZ 
CRIP(2,1) =CZ21—-(CDOIT*€bp21 *C2Z22-CDLI C2 2*C 421 -Cb7 1 Cea 
+CD2i*CD22*CZ11) *CRDGRZ 
Gre(2, 2 )=c422—(CD12 *CO21*CZ22-Cbil 2*CD22° C221 -Cb2z i epez eae 
FeD2 2 COZ2Z C211) CRD Zz 
CRI (1, 3)=(CD12*CZUa-ebiise7i2) *CRErzZ 
CRICl, 4) =(CD1l2*CZ21-CD1LI*C2aZZ een DTZ 
CRIG(2 3) =4ED2 2 *CZ.1 1 =-CD215C7ae Ae eCRDTZ 
CRI? 4 i= (CD22 C221 -Cb21 C222), "CRDTZ 
GRIN3) 1 )=(ebld “221-6 02) C711). *CRDTZ 
GRe(3;, 21 =(eDl2*CZ21-CDz22* CA) ERDEA 
Gk 4., Ly=(eD11 *C222-Cb21*C212) eRpre 
CRI(4, 2) =(CD12*CZ22-Gb22>E212) *CRDTZ 
Shits, 5S) =CZb wen D RZ 
CRI, 4) =€221 *CROTZ 
GRib(4 73) =€7027 CRDTZ 
CR1 (4, 4) =CZ22*CRDTZ 


DO 1300 IQE=0,MXQEG 
IQE1=IQE+1 

£O=2* 105 

De=16 

DQ1=I0+1 
IQX1=4*IO+1 
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TOxZ=TOR Ie 
LOxs=l0xZ2e1 
IQOX4=I0X3 +1 
DO 1100 IPE=0,MXQEG 
IPE1L=IPE+1 
TP=2* TPE 
IPX1=4* IP+1 
IPX2=IPX1+1 
IPX3=IPX2+1 
IPX4=IPX3+1 
DP=IP 
DP1=IP+1 
DBPHI=PISQ* (DP+DQ1) * (DP1-DQ) 
BPHI=4 .D0/DBPHI 
DBZ=PISQ* (DP+DQ1) * (DP1-DQ) * (DP1+DQ1+ONE) * (DP-DQ1) 
BZ=-8 .D0*DQ1/ DBZ 
CRPQ1 (IPX1,1IQX1)=CR1(1,1) *BPHI 
CRPQ1 (IPX2, IOX1) =CR1 (2,1) *BPHI 
CRPQ1 (IPX3,IQX1)=CR1(3,1)*BPHI 
CRPQ1 (IPX4,IQX1)=CR1(4,1)*BPHI 
CREO] (IPX1 , I1OX2) =CR1 (1,2) *BZ 
CRPQ1 (IPX2,IQX2)=CR1(2,2)*BZ 
CRPO1 (IPX3 ,IQX2)=CR1(3,2)*BZ 
CRPOL(IPX4, 1OX2)=CR1(4,2)*BZ 
CRPO1 (IPX1, IOX3) =CR1 (1,3) *BPHI 
CRPOQ1 (IPX2, IQX3) =CR1(2,3)*BPHI 
CRPOQ1 (IPX3, IQX3) =CR1 (3,3) *BPHI 
CRPQ1 (IPX4,1IQX3) =CR1(4,3) *BPHI 
CRPO1 (IPX1, IOX4) =CR1 (1,4) *BZ 
CRPO1 (IPX2,IQX4) =CR1(2,4)*BZ 
CRPO1 (IPX3,1Q0X4) =CR1(3,4)*BZ 
CRPQ1 (IPX4,IQX4) =CR1 (4,4) *BZ 
1100 CONTINUE 
1300 CONTINUE 
DO 2300 I00=0,MXQ0G 
IQ01=1I00+1 
IQ=2*IQ00+1 
IQX1=4*I0+1 
ITOX2=IQX1+1 
EOXRS =10x%2+1 
IOX4=IQ0X3+1 
DQ=I0 
DQO1=1I0+1 
DO 2200 IPO=0,MXQ0G 
IPO1L=IPO+1 
IP=2* IPO+1 
IPX1=4*IP+1 
IPX2=IPX1+1 
IPX3=IPX2+1 
IPX4=IPX3+1 
DP=IP 
DP1=IP+1 
DBPHI=PISQ* (DP+DQ1) * (DP1-DQ) 
BPHI=4.D0/DBPHI 
DBZ=PISQ* (DP+DQ1) * (DP1-DQ) * (DP1+DQ1+ONE) * (DP-DQ1) 
BZ=-8 .DO0*DQ1/DBZ 
GREOI( EPA, LOX) =CR1(1,1)*BPHI 
CRE@ICIPxA2, LOX1)=CR1 (2,1) *BPHI 
GRPOU (GPAs, LOX! )=CR1 (3,1) *BPHI 
CRPO1 (IPX4, IQX1) =CR1(4,1) *BPHI 
CREOMCIPX1 /-TOXx2)=CR1(1,2)*BZ 
CRPOL(IPX2, TOX2) =CR1(2,2)*BZ 
CREOL(CIPXS, [OX2Z) =CR1 (3,2) *BZ 
CRPQ1 (IPX4, IQX2) =CR1 (4,2) *BZ 
CREO PUIEX1, [OX3 )}=CR1 (1,3) *BPHI 
CREGIVLEEX2Z -TOX3) =CR1 (2,3) *BPHI 
CRPQ1 (IPX3, IO0X3) =CR1 (3,3) *BPHI 
CRPQ1 (IPX4,1IQX3) =CR1 (4,3) *BPHI 


iJ) 


CReOlmL(Lex) 7 LOx4) =CR1 (1,4) *BZ 
CRPQ1 (IPX2, IQX4)=CR1 (2,4) *BZ 
CRPQ1 (IPX3, IOX4) =CR1 (3,4) *BZ 
CRPQ1 (IPX4, I1QX4)=CR1 (4,4) *BZ 
2200 CONTINUE 
2300 CONTINUE 
C Set up the matrix Rpq utilized when Z or Delta is nor diagonal, andn < 0. 
CR2 (a), 2) =-Cr1 (1, 2) 
Ck2(2 41) =-Cri(Z, 1) 
CR Oly, )=-CRiIE C1, 3)) 
CR2Z(2774)=—CR1 (2,4) 
CR2(3,1)=-CR1 (3,1) 
CR2Z (4,2) =-CR1L (4, 2) 
CR2 (3,4)=-CR1 (3,4) 
CR2 (4,3) =-CR1 (4,3) 
CRZel ee =CRL (1,1) 
CRZ(2).2) -CRL (2,52) 
CR2 (1,4) =CR1 (1,4) 
CR21@2 55) =GR (2.3) 
CR2 (32,2) =CR1 (3, 2) 
CR2 (4,1)=CR1 (4,1) 
CRZ.055,3 V=CRI(3, 3) 
CR2 (4,4)=CR1 (4, 4) 
DO 3300 IQE=0,MXQEG 
IQE1=IQE+1 
I190=2710E 
DQ=I0 
DQ1=IQO+1 
TQX1=4*IQ+1 
IQX2=IQX1+1 
IQX3=IQX2+1 
IQOX4=IQX3+1 
DO 3100 IPE=0,MXQEG 
IPE1=IPE+1 
IP=2* IPE 
IPX1=4*IP+1 
IPX2=IPX1+1 
IPX3=IPX2+1 
I PX4=IPX3+1 
DP=IP 
DP1=IP+1 
DBPHI=PISQ* (DP+DQ1) * (DP1-DQ) 
BPHI=4.D0/DBPHI 
DBZ=PISQ* (DP+DQ1) * (DP1-DQ) * (DP1+DQ1+ONE) * (DP-DQ1) 
BZ=-8.D0*DQ1/DBZ 
CRPQ2 (IPX1,IQX1)=CR2(1,1) *BPHI 
CRPQ2 (IPX2 ,1QX1) =CR2 (2,1) *BPHI 
CRPQ2 (IPX3 , IQX1) =CR2 (3,1) *BPHI 
CRPQ2 (IPX4,1IQX1)=CR2 (4,1) *BPHI 
CRPOQ2 (IPX1,1QX2) =CR2 (1,2) *BZ 
CGRPO2 (I PX2 -TOXZ)=CR2Z (2,2) “BZ 
CRPQ2 (IPX3 , IQX2) =CR2 (3,2) *BZ 
CRPQ2 (IPX4, IQX2) =CR2 (4, 2) * BZ 
CRPOQ2 (IPX1, IQX3)=CR2 (1,3) *BPHI 
CRPQ2 (IPX2, IQX3) =CR2 (2,3) *BPHI 
CRPQ2 (IPX3 , IOX3) =CR2 (3,3) *BPHI 
CRPQ2 (IPX4, IOX3) =CR2 (4,3) *BPHI 
CRPO2 (IPX1,1IQX4)=CR2 (1, 4) *BZ 
CRPQ2 (IPX2, IOX4) =CR2 (2,4) *BZ 
CRPQ2 (IPX3 , I10X4) =CR2 (3,4) * BZ 
CRPQ2 (IPX4, IQX4)=CR2 (4,4) * BZ 
3100 CONTINUE 
3300 CONTINUE 
DO 4300 IQ0=0,MXQO0G 
T1Q01=I00+1 
TQ=2*IO0+1 
TOX1=4* IQO+1 
IQX2=IQOX1+1 
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LOXxS=TOxZ +E 
TOXx4=10OxX3 +1 
DQ=IQ 
DQ1=I10+l1 
DO 4200 IPO=0,MXQO0G 
IPO1=IPO+1 
IP=2* IPO+1 
IPX1=4*IP+1 
IPX2=IPX1+1 
IPX3=IPX2+1 
IPX4=IPX3+1 
DP=IP 
DP1=IP+1 
DBPHI=PISQ* (DP+DQ1) * (DP1-DQ) 
BPHI=4/DBPHI 
DBZ=PISQ* (DP+DQ1) * (DP1-DQ) * (DP1+DQ1+ONE) * (DP-DQ1) 
BZ=-8.*DQO1/DBZ 
CRPQ2 (IPX1,1IQX1) =CR2 (1,1) *BPHI 
CRPO2 (IPX2, IQXK1) =CR2 (2,1) *BPHI 
CRPQ2 (IPX3 , IQX1) =CR2 (3,1) *BPHI 
CRPQ2 (IPX4, I190X1) =CR2 (4,1) *BPHI 
CRPQ2 (IPX1, IQX2) =CR2 (1,2) *BZ 
CRPQ2 (IPX2, IQX2) =CR2 (2,2) *BZ 
CRPO2 (I PXS , 1OX2Z) =CR2 (3,2) *BZ 
CRPQ2 (IPX4,1QX2) =CR2 (4,2) *BZ 
CRPQ2 (IPX1, 1I0X3) =CR2 (1,3) *BPHI 
CRPQ2 (IPX2,1QX3) =CR2 (2,3) *BPHI 
CREO? (1PxX3, TOXS ) =CR2(3,3) *BPHI 
CRPQ2 (IPX4,10QX3) =CR2 (4,3) *BPHI 
CRPQ2 (IPX1, IQ0X4) =CR2 (1,4) * BZ 
CRPOZ(TPX2, 1O%4)=CR2 (2,4) *BZ 
CRPQ2 (IPX3, 10X4) =CR2 (3,4) *BZ 
CRPQ2 (IPX4, 10X4) =CR2 (4,4) *BZ 
4200 CONTINUE 
4300 CONTINUE 
C Prepare a mtarix for computing the scattering currents on the outer 
C and the inner surfaces 
Pe eclho oNeE. LL) THEN 
RETURN 
END IF 
CRS(1 7.1) =—-CRi(4, 1) 
Gre Gl, Z2j==-CRI (4, 2) 
CR3 (1,3) =-CR1 (4,3) 
CRS Gia) =—-CRiaad 4) 
Che (5s , Doeki(2, 1) 
CRS (3372) =CRI (2.,.2) 
CR3(3,3)=CR1 (2,3) 
CR3(3,4)=CR1 (2,4) 
Creo) =CRL (3,1) 
GRE (2,2)-CR1 (372) 
Ghat 2, 5) =Ch i555) 
CRs (274)-=CR1 (3,4) 
Gre (4, 1)=-CR1 (1,1) 
CR3 (4,2) =-CR1 (1, 2) 
Gro 4,5)=-CR1 (1,3) 
CR3 (4,4) =-CR1 (1,4) 


DO 5200 IQE=0,MXQEG 
IQE1L=IQE+1 
IQ=2* IQE 
DQ=1I9 
DQ1=I0+1 
TOX1=4*I0+1 
IQX2=IQOX1+1 
IQOX3 =IQOX2+1 
IQX4=IOX3+1 

DO 5100 IPE=0,MXQEG 
IPEL=IPE+1 


jz 


1G) 5 — PANG) sho 

IPX1=4*IP+1 

IPX2=I PX1+1 

IPX3=IPX2+1 

IPX4=IPX3+1 

DP=LF 

DPI—1 P+ 

DePHt=FPISo~ (Pr +bOl) * (DP1-DO) 
BPHI=4./DBPHI 

DBZ=PISQ* (DP+DQ1) * (DP1-DQ) * (DP1+DQ1+ONE) * (DP-DQ1) 
BZ=-8.*DQ1/DB2 

CREOo (Lex, TOX)) =(CONE+CR3 (1, 1) )“BPHI 
CRPO31 (IPX2, IQX1) =CR3 (2,1) *BPHI 
CREOOP (UP x3S  @Onl) =CR3 (3,41) “BEAL 

CRPQ31 (IPX4, IOX1) =CR3 (4,1) *BPHI 

GREOS) (TPA) 1042) =CR3 (1,2) ~BzZ 

CRPQ31 (IPX2, IQX2) = (CONE+CR3 (2, 2) ) *BZ 
GREOs TM (TPxs LOX2Z) =CR3 (3,2) *B2 

CRPOS1 (1 PX4, 10X2) =CR3 (4, 2) *BZ 

CREGS 1 (tex, 2OXx3) =CR3 (1,35) “BPH 
CREOSIT(2PAZ, LOX3 ) =CR3 (2,3) “BPHL 

CREOS 1] (1PA37,10X%3) = (CONE+CRS (3.75)))) *BPHI 
CRPO31(1PX4, 10X3) =CR3 (4,3) *BPHI 

CREOS 1] (tPX1, 1OX4) =CR3 (1,4) *Bz 

CRPQ31 (IPX2, IOX4) =CR3 (2,4) *BZ 

CREGS 2 (2X3, [OX4) =CRS (3,4) *B2 

CRPQ31 (IPX4, IQX4) = (CONE+CR3 (4,4) ) *BZ 


GReos2 (2PX1, 1OX1) = (CONE-CR3(171)) *BPHI 

CREOS2 (1 PXK2, [OX1) =-CR3 (2, 1) *BPHI 

CREO i2 EP x3 , 1OX1) =-CR3S (3,1) *BPHL 

CRPOQO32 (IPX4, IQX1) =-CR3 (4,1) *BPHI 

CREOsZa wrx TOK2) =-CR3 (1, 2) 752 

GREQS2 (1PX2, [0OK2) —(CONE-CR3 12, 2) )*BZ 

CREOs2 (TL PX3 , 1OX2)=-CR3(3 42) 2B2 

GRRGs2Z(1PX4, 1OX2Z)=-CR3 (252) *BZ 

GREOs2 rx, 1OX3) =-CR3 (1, 3) “BPHI 

CRPQO32(IPX2,IQX3)=-CR3 (2,3) *BPHI 

GREOS 2 (LT PX37 1OXs ) = (CONE-CRS (37 3) ) *SPHI 

GRPO32 (1 PX4, 1OX3) =-CR3 (4,3) “BPHE 

CRPQ32 (IPX1,IQX4) =-CR3 (1,4) *BZ 

CRPQ3 2 (IPX2, IOX4) =-CR3 (2, 4) *BZ 

CReOS2 (2PXS, 10X44) =-CRS (3,4) "BZ 

CRPQ32 (IPX4,1QX4) =(CONE-CR3 (4,4))*BZ 
CONTINUE 

CONTINUE 


DO 6200 IQ0=0,MXQ0G 
IQ01=IQ00+1 
IQ=2* IQO+1 
DQO=I0 
DQO1=I10+1 
IQX1=4*I0+1 
IQOX2=IQX1+1 
IQOX3=IQX2+1 
IQOX4=IOX3+1 

DO 6100 IPO=0,MXQO0G 
IPE1=IPO+1 
IP=2*ITPO+1 
IPX1=4*IP+1 
IPX2=IPX1+1 
IPX3=IPX2+1 
IPX4=IPX3+1 
DP=IP 
DP1=IP+1 
DBPHI=PISQ* (DP+D0Q1) * (DP1-DQ) 
BPHI=4./DBPHI 
DBZ=PISQ* (DP+DQ1) * (DP1-DQ) * (DP1+DQ1+ONE) * (DP-DQ1) 
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BZ=-8 . *DO1/DBZ 

CRPQ31 (IPX1,IQX1) = (CONE+CR3 (1,1) )*BPHI 
CRPQ31 (IPX2,1IQX1)=CR3 (2,1) *BPHI 

CRPQ31 (IPX3, IQX1) =CR3 (3,1) *BPHI 

CRPQ31 (IPX4,1IQX1) =CR3 (4,1) *BPHI 

CRPQ31 (IPX1, 1IQX2)=CR3 (1,2) *BZ 

CRPQ31 (IPX2, IQX2) = (CONE+CR3 (2,2))*BZ 
CRPQ31 (IPX3, 1QX2) =CR3 (3,2) *BZ 

CRPQ31 (IPX4,1IQX2)=CR3 (4,2) *BZ 

CRPQ31 (IPX1,1QX3)=CR3 (1,3) *BPHI 

CRPQ31 (IPX2,1IQX3)=CR3 (2,3) *BPHI 

CREO’ TPX3, 1Ox%3) = (CONE+CR3 (3,3) ) *BPHI 
CRPQ31 (IPX4,1I0QX3) =CR3 (4,3) *BPHI 

CRPQ31 (IPX1,IQX4)=CR3 (1,4) *BZ 
CRPQ31(IPX2,10X4) =CR3 (2,4) *BZ 

CRPQ31 (IPX3,1QX4) =CR3 (3,4) *BZ 
CRPQ31(IPX4,1IQX4) =(CONE+CR3 (4,4)) *BZ 


cS 
CRPQ32 (IPX1, IQX1) = (CONE-CR3 (1,1) ) *BPHI 
CRPQ32 (IPX2,IQX1)=-CR3 (2,1) *BPHI 
CRPQ32 (IPX3, IQX1) =-CR3 (3,1) *BPHI 
CRPOQ32 (IPX4,IQOX1) =-CR3 (4,1) *BPHI 
CREO¢ZtLex LL. £022) =—-CR3 (1,2) *BZ 
CREQ22) ERX? , £042) =(CONE-CR3 (2,2))*BZ 
CRPQ32 (IPX3 , IQX2) =-CR3 (3, 2) *BZ 
CRPQ32 (IPX4,1IQX2)=-CR3 (4,2) *BZ 
CRPQ32 (IPX1, IQX3) =-CR3 (1,3) *BPHI 
CRPQ32 (IPX2, IQX3) =-CR3 (2,3) *BPHI 
CREQS7Z (lr xs LOX3 ) =(CONE-CR3 (3, 3) )* BPHE 
CRPQ32 (IPX4, I9X3) =-CR3 (4,3) *BPHI 
CRPOQ32 (IPX1,IQX4)=-CR3 (1,4) *BZ 
CRPQ32 (IPX2 , I9X4) =-CR3 (2,4) *BZ 
CRPOS2 (1 PX3, LOx4) =—Grei3 , 4) *BZ 
CRPQ32 (IPX4, I1QX4) = (CONE-CR3 (4,4) ) *BZ 

6100 CONTINUE 

6200 CONTINUE 
RETURN 

9001 FORMAT('The given inside and outside impedance have a singular', 
+ ' sum. Execution is terminated. ') 
END 


H. SUBROUTINE GNPQFN 


SUBROUTINE GNPOFN 
CREEK RK KKK KK KK RK RK KKK KEK KR KKK KKK KKK KKK KKK KKK RK KERR KEKE KK KEKE KKK KKK KR KKK KK KKK KKK 
INCLUDE 'REALTP.INC' 
INCLUDE '‘MAINDM.INC' 
COMMON /GCONST/ DKH, DKA, HH,HA,HSQ, DASO, HSQN, DASQN, HHSQ, HDASQ, 
+ RAH, RAHSQ, DHA, DAH 
CHARACTER FILEEVEN1*12, FILEODD1*12 
C Set up the file names 
NC1=INT (DKA) 
DE=NCr 
NC=100000.*DKH+NC1 
DC1=DKA-DC 
NC2=INT (1000.*DC1) 
C Set up the indices of file names of Green's function, and check whether 
C these files exist. If these files exist, then use it directly. Otherwise, 
C call subroutine XPQINI1. 
FILEEVENIl='E : : 
FILEODD1='0O : 
IGE=8*2* (MAXPEG+1) * (MAXPEG+2) 
IGO=8*2* (MAXPOG+1) * (MAXPOG+2) 
WRITE “(FILEEVENIA2:5), "(1/7 -7)*) Ne¢ 
WRITE (ELILEODDI (248), (17.7) ") NE 
WRITE, {FILEEVENL(TO-12),  (i3.3)" ) Nez 
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Veeepeeer LEEODD] (10212), * (13.3) '*) NC2 
OPEN (28,ACCESS='DIRECT' , FILE=FILEEVEN1, RECL=IGE, IOSTAT=I0S, 


= STATUS='OLD' ) 
feowtoOs -NE. 0) THEN 
CLOSE (28) 
CALL XPQINI1 (DKH, DKA) 
ELSE 
OPEN (29,ACCESS='DIRECT' , FILE=FILEODD1, RECL=IGO, IOSTAT=I10S, 
+ STATUS='OLD' ) 
TF) (10S. -NE2 0) THEN 
CEOSE. (29) 
CALL XPQINI1 (DKH, DKA) 
ELSE 
CLOSE (28) 
CEOSE (29) 
CALL XPQINI (DKH, DKA) 
END IF 
ENDIF 
RETURN 
END 


Cc 
SUBROUTINE XPQINI (DKHIN, DKAIN) 
CREEK KEKKKEKKKRK EKER KK RHEE EK KREKRKREKKEKK KEK KK KKK KEKE KKKKKEKEKKEKKEKK EK KKK KK KKK 
TINGLUDE “REALTP:. INC” 
ENGLUDE  “CMEXTP INC” 
INCLUDE 'MAINDM.INC' 
DIMENSION GNE(4, (MAXPEG+1) * (MAXPEG+2)/2), 


+ GNO (4, (MAXPOG+1) * (MAXPOG+2) /2) 
DIMENSION CGNE(0:MAXPEG, 0:MAXPEG,KNDIM1+1), 
= CGNO (0:MAXPOG, 0: MAXPOG, KNDIM1+1) 
DIMENSION CDGNE(0:MAXPEG, 0:MAXPEG, KNDIM1+1), 
+ CDGNO (0:MAXPOG,0:MAXPOG, KNDIM1+1) 
COMMON /CRNTDM/ NMAX, MXNG, IOMAX, IQMAX1, IOQMAX2, IXCRNT, ICRNT,MXQEG, 
+ MXOO0G 


COMMON /GPOQTMP/ CGNE,CDGNE, CGNO, CDGNO 
SAVE /GPQTMP/ 


CHARACTER FILEEVEN1*12, FILEODD1*12 
DKH=DKHIN 
DKA=DKAIN 
S 
C Initialize the matrix CGNE, CGNO, CDGNE, CDGNO 
DOV Sa 1C— 71 KNDIMI+1 
DO 102 IB= 0,MAXPEG 
DO 101 IA= 0,MXPEG 
CGNE(IA,IB, IC) =CZERO 
CDGNE (IA, IB, IC) =CZERO 
slows CONTINUE 
102 CONTINUE 
DO 104 ID=0,MAXPOG 
DO 103 IE=0,MAXPOG 
CGNO(IE, ID, IC) =CZERO 
CDGNO (IE, ID, IC) =CZERO 


103 CONTINUE 
104 CONTINUE 
105 CONTINUE 

C 


C Set up the file name 
N@l=INT (DKA) 
DE=NCI 
NC=100000.*DKH+NC1 
DC1=DKA-DC 
NC2=INT(1000.*DC1) 
FILEEVENI1='E : ; 
FILEODD1='0 : 
IGE=8*2* (MAXPEG+1) * (MAXPEG+2) 
IGO=8*2* (MAXPOG+1) * (MAXPOG+2) 
WRITE (FILEEVENI (2:8), ° (1727).') NC 
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300 


400 
500 


600 


700 
800 
900 


as 


WRITE (FILEODD] (2:8 )jaGl7. 7)? jae 

WRITE (PILEEVEND (10:12) Si 3e3). tee 

WRITE SUFI LEODDI (10. 12). Geseeo) Ne2 

OPEN (28,ACCESS='DIRECT' , FILE=FILEEVEN1, RECL=IGE, IOSTAT=I0OS, 
STATUS='OLD' ) 


OPEN (29,ACCESS='DIRECT' , FILE=FILEODD1, RECL=IGO, IOSTAT=IOS, 
STATUS='OLD' ) 


DO 900 NI=1,MXNG+1 
The following values N, DN and DNH are passed to the G-computation 
related subroutines through the common block /NCONST/. 
READ (28,REC=NI) GNE 
IRECE=0 
DO 500 IQE=0,MXQEG 
DO 300 IPE=0,IQE-1 
CGNE (IPE, IQE, NI) =CGNE (IQE, IPE,NI) 
CDGNE (IPE, IOQE,NI)=CDGNE(IQE, IPE,NI) 
CONTINUE 
DO 400 IPE=IQE,MXQEG 
IRECE=IRECE+1 
GR=GNE (1, IRECE) 
GI=GNE (2, IRECE) 
GDR=GNE (3, IRECE) 
GDI=GNE (4, IRECE) 
CGNE (IPE, IQE,NI) =DCMPLX (GR, GI) 
CDGNE (IPE, IQE, NI) =DCMPLX (GDR, GDI) 
CONTINUE 
CONTINUE 
READ (29,REC=NI) GNO 
IRECO=0 
DO 800 IQ0=0,MXQ0G 
DO 600 IPO=0,1900-1 
CGNO (IPO, IQ0, NI) =CGNO(I0Q0, IPO, NTI) 
CDGNO (IPO, IQ0,NI) =CDGNO(IQ0, IPO,NTI) 
CONTINUE 
DO 700 IPO=IQ0,MXQ0G 
IRECO=IRECO+1 
GR=GNO (1, IRECO) 
GI=GNO (2, IRECO) 
GDR=GNO (3, IRECO) 
GDI=GNO (4, IRECO) 
CGNO (IPO,1IQ0,NI) =DCMPLX(GR,GI) 
CDGNO (IPO, IQ00, NI) =DCMPLX (GDR, GDI) 
CONTINUE 
CONTINUE 
CONTINUE 
CLOSE (23) 
CLOSE (29) 
RETURN 
END 


SUBROUTINE XPQINI1 (DKHIN, DKAIN) 


CARER KK KKK KKK kk KK KKkk KK KKK Kk kkk Kk Kk Keke Kk KKK KKK KKK KKK KKK KKK KKK KKK KKK 


+ 


ay 


+} 


ot 


INCLUDE ‘REALTP.INC' 

INCLUDE ‘CMPXTP.INC' 

INCLUDE 'MAINDM.INC' 

INCLUDE ‘GPQNDM.INC' 

INGEUDES LEMITS.INC * 

DIMENSION GNE (4, (MAXPEG+1) * (MAXPEG+2)/2), 
GNO (4, (MAXPOG+1) * (MAXPOG+2) /2) 

DIMENSION CGNE (0:MAXPEG, 0:MAXPEG, KNDIM1+1), 
CGNO (0:MAXPOG, 0: MAXPOG, KNDIM1+1) 

DIMENSION CDGNE(0:MAXPEG,0:MAXPEG,KNDIM1+1), 
CDGNO (0: MAXPOG, 0: MAXPOG, KNDIM1+1) 

COMMON /CRNTDM/ NMAX,MXNG, IQMAX, IOQMAX1, IOQMAX2, IXCRNT, ICRNT, MXQEG, 

MXQOG 
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COMMON /GCONST/ DKH, DKA,HH,HA,HSQ, DASQ, HSQN, DASQON, HHSQ, HDASQ, 
+ RAH, RAHSQ, DHA, DAH 

COMMON /SUMLMT/ DHMAX, DAMAX, DSNG, IHMAX, IAMAX, ISNG 

COMMON /NCONST/ DN, DNH,N 

COMMON /GPOTMP/ CGNE,CDGNE,CGNO,CDGNO 

GHARACTER FILEBVENI*12, FILEODDI*12 

SAVE /GPOQTMP/,/GCONST/, /SUMLMT/ 


DKH=DKHIN 
DKA=DKAIN 
C Initialize the matrix CGNE, CGNO, CDGNE, CDGNO 
HOBLOS 1¢G= 1, KNDEMTs1 
DO 102 IB= 0,MAXPEG 
DOs Ol TA=305 MXPEG 
CGNE (IA, IB, IC) =CZERO 
CDGNE (IA, 1B, 1€)=CZERO 
Oa CONTINUE 
LO2 CONTINUE 
DO 104 ID=0,MAXPOG 
DO 103 IE=0,MAXPOG 
GCGNO(IE, ED, 1C) =CZERO 
CDGNO (IE, ID, IC) =CZERO 


103 CONTINUE 
104 CONTINUE 
iOS CONTINUE 

c 


C Determine the maximum number of terms for r- and m- series sums and 
C pass them through /SUMLMT/: 
3 


REFC=FZERO* LOG (TWO) -LOG (PI2) -ONE 
REFH=HALF* (REFC-LOG (DKH) ) 
REFA=HALF* (REFC-LOG (DKA) ) 


IQMX2=IQMAX+2 
DHMAX=IQMX2 
DO 1100 WHILE ( (DHMAX+HALF) * (LOG (DHMAX/DKH)+ONEN) .LT. REFH) 
DHMAX=DHMAX+ONE 
1100 CONTINUE 
IHMAX=INT (DHMAX) 


& 
DAMAX=AINT (DKA) +ONE 
DO 1200 WHILE ( (DAMAX+HALF) * (LOG (DAMAX/DKA)+ONEN) .LT. REFA) 
DAMAX=DAMAX+ ONE 
1200 CONTINUE 
C 
IAMAX=INT (DAMAX) 
Cc 


ISNG=INT (TWO* DKH+ONE-EPS8) 
ISNG=MAX (IFPBIT, ISNG, IOMX2) 
C Checking dimensions. 
IF (IHMAX .GT. MXMREG) THEN 
WRITE (*,*) ‘Warning: IHMAX = ', IHMAX,' > MXMREG = ',MXMREG 
WRITE(*,*) '‘IHMAX IS SET TO MXMREG IN XPOINI1' 
IHMAX=MXMREG 
END IF 
ISN1=MXMSNG-N-1 
EP{(ISNG .GT. ISN1) THEN 


WRITE(*,*) ‘Warning: ISNG = ',ISNG,' > MXMSNG-N-1 = ',ISN1 
WRITE(*,*) ‘ISNG IS REDUCED TO MXMSNG-N-1 IN XPQINI1’ 
ISNG=ISN1 
HA) OO Oy 
(& 
DSNG=ISNG 
Cc 


FILEEVEN1='E ‘ ’ 
FILEODD1='0O 5 ; 
NC1=INT (DKA) 

DC=NC1 
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NC=DKH*100000.+NC1 
DC1=DKA-DC 

NCZ=INT (1000. * Del) 

IGE=8*2* (MAXPEG+1) * (MAXPEG+2) 
IGO=8*2* (MAXPOG+1) * (MAXPOG+2) 


WRITE (FILEEVENI( 228) 51727) )) Ne 
WRITE (FILEODDI(228)2 (19.2) ) NC 
WRITE. (FILEEVEN] (0212), (13.3) ") NeZ2 
WRITE (FE TLEODD we ye (13.3) ) NC2 


OPEN (28,ACCESS='‘ DIRECT’ , FILE=FILEEVEN1 , RECL=IGE, STATUS='UNKNOWN' ) 

OPEN (29,ACCESS='DIRECT' , FILE=FILEODD1, RECL=IGO, STATUS=' UNKNOWN ' ) 
C Initialize CGNE, CDGNE, CGNO and CDGNO 

DO 1900 NI=1,MXNG+1 


e The following values N, DN and DNH are passed to the G-computation 
Cc related subroutines through the common block /NCONST/. 

N=NI-1 

DN=N 

DNH=DN+HALF 

IRECE=0 

DO 1500 IQE=0,MXQEG 
IQ=2* IQE 
DO 1300, TPE=0) TOE-1 
IP=2* IPE 


CGNE (IPE, IQE,NI) =CGNE(IQE, IPE,NI) 
CDGNE (IPE, IQE, NI) =CDGNE (IQE, IPE, NI) 
13/00 CONTINUE 
DO 1400 IPE=IQE,MXQEG 
IRECE=IRECE+1 
IP=2*IPE 
CALL GDN(IP,1I0,GI,GDI,GR, GDR) 
GNE (1, IRECE) =GR 
GNE (2, IRECE) =GI 
GNE (3, IRECE) =GDR 
GNE (4, IRECE) =GDI 
CGNE (IPE, IQE, NI) =DCMPLX (GR, GI) a 
CDGNE (IPE, IOQE, NI) =DCMPLX (GDR, GDI) 
1400 CONTINUE 
500 CONTINUE 
WRITE (28,REC=NI) GNE 
IRECO=0 
DO 1800 IQ0=0,MX00G 
IQ=2*IQO0+1 
DO 1600 IPO=0,1I00-1 
IP=2*IPO+1 
CGNO (IPO, I1Q0, NI) =CGNO(IQ0, IPO,NI) 
CDGNO (IPO, 1IQ0,NI) =CDGNO(IQO, IPO,NI) 
1600 CONTINUE 
DO 1700 IPO=IQ0, MXQOG 
IRECO=IRECO+1 
TP=2* IPO+1 
CALL GDN(IP,10,GI,GDI,GR, GDR) 
GNO (1, IRECO) =GR 
GNO(2,IRECO)=GI 
GNO (3, IRECO) =GDR 
GNO (4, IRECO) =GDI 
CGNO (IPO, IQ0,NI) =DCMPLX (GR, GI) 
CDGNO (IPO, I00, NI) =DCMPLX (GDR, GDI) 
1700 CONTINUE 
1800 CONTINUE 
WRITE (29,REC=NI) GNO 
1900 CONTINUE 
CLOSE (28) 
CLOSE (29) 
RETURN 
END 


SUBROUTINE GDN(IPIN, IQIN, GIOUT, GDIOUT, GROUT, GDROUT) 


CEKRREKREREKKEKEKEEREKEKREKRERKREERKEEKRAEERKRRRERREEAARRKRERRER YARRA A A 
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C This subrouitne sets up P and Q dependent constants and passes along 
C /PCONST/, then calls subroutines GDREG and GDSNG to compute G(P,0,N) 
C and its derivative. 
CREEK KKK KEKE RK KKK KKK KEKE KEKKE KK KEKE KK KKK KEKE KKK KEK KEKE KKK EK KEKE KK KEKKEK KEKE K 
e 
INCLUDE 'REALTP.INC' 
COMMON /PCONST/ DP,DQ,DS,DD,DSH,DDH,DSSQ,DDSO,IP,IQ,IS,ID 
C Check and transform input variables. 
IP=IPIN 
IQ=IQIN 
ISPQ=IP+IQ 
IS=ISPQ/2 
IF ((ISPQ+1)/2 .GT. IS) THEN 
GIOUT=ZERO 
GDIOUT=ZERO 
GROUT=ZERO 
GDROUT=ZERO 
WRITE Ge 9001). £Seo 
9001 FORMAT('The parameters P and Q have a sum of the ODD integer ',I4, 
+/,', G(P, 9, N) and its derivative have been set to 0.°) 
RETURN 
END? iF 
IF (IP .LT. IQ) THEN 
IP=IQOIN 
IQ=IPIN 
END IF 


DP=IP 
DQ=1IQ 


ID=IS-=10 
DS=IS 

DD=ID 
DSH=DS+HALF 
DDH=DD+HALF 
DSSQ=DS*DS 
DDSQ=DD* DD 


CALL GDREG(GI,GDI,GRR,GDRR) 

CALL GDSNG (GRS,GDRS) 

GIOUT=GI 

GDIOUT=GDI 

GROUT=GRR+GRS 

GDROUT=GDRR+GDRS 

RETURN 

END 
CC 

SUBROUTINE GDREG(GIOUT, GDIOUT, GROUT, GDROUT) 
CEEKKKREEKKRKKKEEKKEKKKKKEKEKKKKEKEKKKEKKEKKKKEKKKKEKEKKEKKKKEKKKEKRKKEKKEKEKEKEKEE 
C This subrouitne computes the regular part of G(P,Q,N) and its 
C derivative. 
CER R KEK REKEEKKKEKKEEKEKKKKEKKEEKKKEREEKEEKEKEEEKKEKKEKEKKEEKRKEKEEKKKRKEKEKEKEEKEKKEKEKE 
% 

INCLUDE 'REALTP.INC' 
INCLUDE 'GPONDM.INC' 

COMMON /GCONST/ DKH,DKA,HH,HA,HSQ, DASQ, HSON, DASON, HHSQ, HDASQ, 

+ RAH, RAHSQ, DHA, DAH 

COMMON /NCONST/ DN,DNH,N 

COMMON /PCONST/ DP,DQO, DS, DD,DSH, DDH, DSSO, DDSO_, IP, 19,1S,20 

COMMON /SUMLMT/ DHMAX, DAMAX, DSNG, IHMAX, IAMAX, ISNG 

SAVE /GCONST/, /SUMLMT/ 


C 

C Reserve working space to store r-independent numbers. 
DIMENSION GIM(MXMREG) , GRRM (MXMREG) 

Cc 

C Computation starts. 

6: 

C Compute overal constant factors: 
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300 


400 


c 


GRRF=QUAR* DKH/ PISQ/ (DSSQ-QUAR) / (QUAR-DDSQ) / (DN+ONE) 
GIF=HALF/DSH 
DIJIN=ZERO 
DO 300 JN=1,N 
DJN=DJN+ONE 
HATN=HA/DJN 
GRRF=GRRF* HATN* HATN 
GIF=GIF*HATN*HA/ (DJN+DSH) 
CONTINUE 
DJP=ZERO 
DO 400 JQ=1,10 
DJP=DJP+ONE 
GIF=(HHSQ/DJP/DJP) *GIF 
CONTINUE 
BOe500 Ur=lO+1,IP 
DJP=DJP+ONE 
GIF=(HH/DJP) *GIF 
CONTINUE 
IF ((ID+1)/2 .GT. ID/2) GIF=-GIF 


C Compute GI, GDI, GRREG, and GDRREG. 


e 


600 


700 


Compute r-independent factors and store in GIM and GRRM: 
SMH=DSH+ONEN 
SM1=DS 
SM2=DS+DS 
PM1=DP 
QM1=D0 
THRH=HALF 
DM1=ZERO 
THRS=HALF+DS 
THRD=HALF+DD 
THRDN=HALF-DD 
THRSN=HALF-DS 
DO 600 JM=1,IHMAX 
SMH=SMH+ONE 
SM1=SM1+ONE 
SM2=SM2+ONE 
PM1=PM1+ONE 
QM1=QM1+ONE 
THRH=THRH+ONE 
DM1=DM1+0ONE 
THRS=THRS+ONE 
THRD=THRD+ONE 
THRDN=THRDN+ONE 
THRSN=THRSN+ONE 
GIM (JM) = (SMH/SM2) * (SMH/PM1) * (SM1/0M1) * (HSQON/DM1) 
GRRM (JM) = (THRH/THRS) * (HSON/THRD) * (DM1/THRDN) * (DM1/THRSN) 
CONTINUE 


Compute r- and m-sum for GI and GRR. 
Setup initial r related values 
DJR=DAMAX 
DNR=DN+DJR 
DNHR=DNR-HALF 
D2NR=DN+DNR 
DNSHR=DNR+DSH 
DN2R=DNR+ONE 
m-sum for r=IAMAX 
DNSHRM=DNSHR+DHMAX 
DN2RM=DN2R+DHMAX 
ANNEXTR=ONE+GIM (IHMAX) /DNSHRM 
ANROR=ONE+GRRM (IHMAX) /DN2RM 
DO 700 JM=IHMAX-1,1,-1 
DNSHRM=DNSHRM+ONEN 
DN2RM=DN2RM+ONEN 
ANNEXTR=ONE+ANNEXTR*GIM(JUM) /DNSHRM 
ANROR=ONE+ANROR*GRRM (JM) /DN2RM 
CONTINUE 


8/ 


ANNEXT=ANNEXTR 
DANNEXT=DNR* ANNEXTR 
ANRO=ANROR 
DANRO=DNR* ANROR 
DO 900 JR=IAMAX,1,-1 
Cc Compute factors for r=JR 
ATMP= (DNHR/D2NR) * (DASON/DJR) 
ATMPI=ATMP/DNSHR 
ATMPR=ATMP/DN2R 
Cc Setup new r related values for r=JR-1 
DJR=DJR+ONEN 
DNR=DNR+ONEN 
DNHR=DNHR+ONEN 
D2NR=D2NR+ONEN 
DNSHR=DNSHR+ONEN 
DN2R=DN2R+ONEN 
Cc m-sum for r=JR-1 
DNSHRM=DNSHR+DHMAX 
DN2RM=DN2R+DHMAX 
ANNEXTR=ONE+GIM ( IHMAX) / DNSHRM 
ANROR=ONE+GRRM (IHMAX) /DN2RM 
DO 800 JM=IHMAX-1,1,-1 
DNSHRM=DNSHRM+ONEN 
DN2RM=DN2 RM+ONEN 
ANNEXTR=ONE+ANNEXTR*GIM(JM) /DNSHRM 
ANROR=ONE+ANROR*GRRM(JM) /DN2RM 
800 CONTINUE 
ANNEXT=ANNEXTR+ANNEXT* ATMPI 
DANNEXT=DNR* ANNEXTR+DANNEXT * ATMPI 
ANRO=ANROR+ANRO* ATMPR 
DANRO=DNR* ANROR+DANRO*ATMPR 
900 CONTINUE 


GIOUT=ANNEXT* GIF 
GDIOUT=DANNEXT*GIF 
GROUT=ANRO*GRRF 
GDROUT=DANRO *GRRF 


RETURN 

END 
c 

SUBROUTINE GDSNG(GROUT,GDROUT) 
CER KKKRAKKKREKREKREEREKKEKKEEKEKREKKKKEKKEKEKRKKEKKKKEKAKRKKKKEKEKKEKKKKKEKKKKEKK K 
C This subrouitne computes the 'Singular' part of G(P,Q,N) and its 
C derivative. 
CER EK REAR REKEKEEKKEEERREKKKKEKEKKEREKKEKEKKKEKEEKKKEKREKEKKEKEEREKEKKEKKEEKEKE 
& 

INCLUDE 'REALTP.INC' 

INCLUDE 'GPONDM.INC' 

COMMON /GCONST/ DKH, DKA, HH, HA,HSQ, DASQ,HSON, DASON, HHSQ, HDASQ, 

+ RAH, RAHSQ, DHA, DAH 

COMMON /NCONST/ DN,DNH,N 

COMMON /PCONST/ DP,DO,DS,DD, DSH; DDH, DSSQO, DDSO, LF, 10, 15, ED 

COMMON /SUMLMT/ DHMAX, DAMAX, DSNG, IHMAX, IAMAX, ISNG 

SAVE /GCONST/, /SUMLMT/ 


C Reserve working space to store r-independent numbers. 
DIMENSION GS1M(MXMSNG) ,GS2M(MXMSNG) , 
+ GR2M(MXMSNG) , GR2R(0:MXMSNG) , GDR2R(0:MXMSNG) 
Bs 
C Computation starts. 
© 


MXGR2M=ISNG+N 


C Compute overal constant factors: 
GRSF=QUAR/ PISQ/DKH 


C Compute GR1, GDR1, GR2, and GDR2: 


88 


Cc Compute initial constants: 


S1LHA=TWO* LOG (HXD/RAH) 
DK=HALF 
S1SD0=ZERO 
DO 1000 JK=1,I1ID 
S1SD0=S1SD0+ONE/DK 
DK=DK+ONE 
1000 CONTINUE 
S1SD0=TWO*S1SD0 
DO 1100 JK=ID+1, IS 
S1SD0=S1SD0+ONE/DK 
DK=DK+ONE 
1100 CONTINUE 
S1SD0=-TWO*S1TSDO0 
SNON=ZERO 
SN10=ZERO 
SN20=ZERO 
IF (N .GT. 0) THEN 
DJN=ONE 
DJNH=HALF 
DO 1200 JN=1,N-1 
DJINI=ONE/DJIN 
DJNHI=ONE/DJNH 
SN10=DJNI* DJINHI+SN10 
SN20=DJNHI*DJNHI* ( (ONE-QUAR*DJNI) *HALF * DJINI+ONE) +SN20 
DJN=DJN+ ONE 
DJNH=DJNH+ONE 
1200 CONTINUE 
DINI=ONE/DJIN 
DJNHI=ONE/DJNH 
SNON=S1LHA+HALF* DJNHI-SN10*QUAR 
SN10=(DJNI* DJINHI+SN10) *QUAR 
SN20=HALF* (DJNHI* DJNHI* ( (QUAR* DJNI+ONEN) *HALF* DJINI+ONEN) -SN20 ) 
END IE 
SN10=S1LHA-SN10 
SN20=PISQ*THIR+SN20 


c Compute r-independent terms and store in GR2M, GS1M, and GS2M: 
GR2JM=ONE 
S1SDM=S1SD0 
S2SDM=ZERO 
DM=ZERO 
DMH=-HALF 
DO 1300 JM=1,MXGR2M 
DM=DM+0ONE 
DMH=DMH+ ONE 
DMI=ONE/ DM 
DMHI=ONE/DMH 
DMHSQ=DMH* DMH 
DMHSS=DMHSQ-DSSQ 
DMHSD=DMHS0Q-DDS0Q 
GR2JM= (DMHSS*DMI* DMI) * (DMHSD*DMI*DMHI) *GR2JM 
GR2M (JM) =GR2JM 
TWSMHS=TWO* DSSQ/DMHSS 
TWDMHD=TWO* DDSQ/DMHSD 
S1SDM=S1SDM- (TWSMHS+TWDMHD+DMI+ONE) * DMHI 
S2SDM=S2SDM- ( (TWSMHS+THR) *TWSMHS+ (TWDMHD+THR) * TWOMHD+ 
+ (ONE-QUAR*DMI) *TWO*DMI+ONE) /DMHSQ 
GS1M (JM) =S1SDM 
GS2M(JM)=S2SDM 
1300 CONTINUE 
Cc 
Ss Compute m-sum for GR1, GDR1 and store in GR2R, GDR2R 
IF (N .GT. 0) THEN 
DJIR=DN+ONEN 
DJRH=DJR-HALF 
SNOR=SNON 


89 


SGR1M=SNOR+S1SD0 
SDGR1M=DJIR* SGR1M+ONEN 
DJM=ONE 
DJIRM=DJIR 
DO 1400 JM=1,N-1 

SR1M=SNOR+GS1M (JM) 
HSQF=HSQ/ DJIM/DIJRM 
SGR1M=SGR1M* HSOQF+SR1M*GR2M (JM) 
SDGR1M=SDGR1M* HSOF+ (DJR* SR1M+ONEN) *GR2M (JM) 
DJIM=DJM+ONE 
DJIRM=DJRM+ONEN 

1400 CONTINUE 
GR2R (N-1) =SGRIM* TWO 
GDR2R (N-1) =SDGR1M* TWO 

DO 1600 JR=N-2,0,-1 
SNOR=SNOR+ONE/DJRH+ONE/ (DN-DJR) -ONE/ (DN+DJR) 
DJR=DJIR+ONEN 
DJRH=DJRH+ONEN 
SGR1IM=SNOR+S1SD0 
SDGR1M=DJIR* SGR1M+ONEN 
DJM=ONE 
DIJRM=DIR 
DO 1500 JM=1,JIR 

SR1M=SNOR+GS1M (JM) 
HSQF=HSQ/DJM/DJRM 
SGR1M=SGR1M* HSQF+SR1M* GR2M (JM) 
SDGR1M=SDGR1M* HSOF+ (DJR* SR1M+ONEN) *GR2M (JM) 
DJIM=DJIM+ONE 
DJRM=DJRM+ONEN 

LSso0 CONTINUE 
GR2R (JR) =SGRIM* TWO 
GDR2R (JR) =SDGR1IM* TWO 


1600 CONTINUE 
END IF 
Cc 
ie Compute m-sum for GR2, GDR2 and store in GR2R, GDR2R 
DRN=DN 
SN1R=SN10 
SN2R=SN20 


SR2MM=SN1R+S1SD0 
SR2RMS=SR2MM* SR2MM+SN2R 
DR2RMS=DRN* SR2RMS-TWO* SR2MM 
DJM=ONE 
DJRM=DRN 
POs 700 JIM=1,N 
SR2MM=SN1R+GS1M (JM) 
SR2RM=SR2MM* SR2MM+SN2R+GS2M (JM) 
DR2RM=DRN* SR2RM-TWO*SR2MM 
HSQF=HSQ/DJM/DJIRM 
SR2RMS=SR2RMS* HSOF+SR2RM*GR2M (JM) 
DR2RMS=DR2RMS* HSQF+DR2RM*GR2M (JM) 
DIM=DJIM+ONE 
DJRM=DJRM+ONEN 
1700 CONTINUE 
GR2R(N) =SR2RMS 
GDR2R(N) =DR2RMS 
DJR=ZERO 
DO 1900 JR=N+1,MXGR2M 
DJR=DJR+ONE 
DRN=DRN+ONE 
DJRI=ONE/DIR 
DRNHI=ONE/ (DRN-HALF) 
DR2NI=ONE/ (DRN+DN) 
DNHRF=DNH* DRNHI* DR2NI 
SN1LR=DJRI-~DNHRF+SN1R 
SN2R=DJRI* DJRI-~- (DRNHI+DR2NI) *DNHRF+SN2R 
SR2MM=SN1R+S1SD0 
SR2RMS=SR2MM* SR2MM+SN2R 
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DR2RMS=DRN*SR2RMS-TWO*SR2MM 
DJM=ONE 
DJRM=DRN 
DO 1800 JM=1,JUR 

SR2MM=SN1R+GS1M (JM) 
SR2RM=SR2MM* SR2MM+SN2R+GS2M (JM) 
DR2RM=DRN* SR2RM-TWO* SR2MM 
HSQF=HSQ/DJM/DJIRM 
SR2RMS=SR2RMS*HSOF+SR2RM*GR2M (JM) 
DR2RMS=DR2RMS*HSOF+DR2RM* GR2M (JM) 
DJM=DJM+ONE 
DJRM=DJRM+ONEN 

1800 CONTINUE 
GR2R(JR) =SR2RMS 
GDR2R (JR) =DR2RMS 


1900 CONTINUE 
C 
c Compute r-sum: 
Cc 
IF (RAHSQ .GT. HALF) THEN 

Cc 
SC Convergence acceleration via Euler-Abel Transfomation: 
Sy 

DJR=ZERO 


DJRH=DJR-HALF 
DJRN=DJR+DN 
DJRNN=DJR-DN 

te Gr. 0) THEN 
RFCTOR=ONE/DN 
GR2R (0) =RFCTOR*GR2R (0) 
GDR2R (0) =RFCTOR*GDR2R (0) 

DOF 2000 JR=1,N-1 
DJR=DJR+ONE 
DJRH=DJRH+ONE 
DJRN=DJRN+ONE 
DJRNN=DJRNN+ONE 
RFCTOR= (DJR* DJRH/DJRN/ DJRNN) *RAHSOQ*RFCTOR 
GR2R (JR) =RFCTOR*GR2R (JR) 
GDR2R (JR) =RFCTOR*GDR2R (JR) 
2000 CONTINUE 

DJR=DJR+ONE 
DJRH=DJRH+ONE 
DJRN=DJRN+ONE 
DJRNN=DJRNN+ONE 
RFCTOR=HALF* (HALF-DN) * RAHSQ* RFCTOR 
GR2R(N) =RFCTOR*GR2R (N) 
GDR2R (N) =RFCTOR*CDR2R (N) 

ELSE 
RFCTOR=ONE 

END IF 

DO 2100 JR=N+1,MXGR2M 
DJR=DJR+ONE 
DJRH=DJRH+ONE 
DJRN=DJRN+ONE 
DJRNN=DJIRNN+ONE 
RFCTOR= (DJR* DJRH/ DJRN/DJRNN) * RAHSQO* RFCTOR 
GR2R (JR) =RFCTOR*GR2R (JR) 
GDR2R (JR) =RFCTOR* GDR2R (JR) 


2100 CONTINUE 
c 
c Compute zeroth-order delta-r term: 


DO 2300 JDELTR=1,MXGR2M 
DO 2200 JR=MXGR2M, JDELTR,-1 
GR2R (JR) =GR2R (JR-1) -GR2R (JR) 
GDR2R (JR) =GDR2R(JR-1) -GDR2R (JR) 
2200 CONTINUE 
2300 CONTINUE 
SGR2R=HALF*GR2R (MXGR2M) 
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SDGR2R=HALF* GDR2R (MXGR2M) 
DO 2400 JR=MXGR2M-1,0,-1 
SGR2R= (SGR2R+GR2R (JR) ) *HALF 
SDGR2R= (SDGR2R+GDR2R (JR) ) *HALF 
2400 CONTINUE 
GROUT=GRSF*SGR2R 
GDROUT=GRSF*SDGR2R 


ELSE 


Direct summation when RAHSQ is no greater than 1/2: 


aaQ0N 


DJR=DSNG 
DJRN=DJR+DN 
DJRNH=DJRN-HALF 
DIRZN=DJRN+DN 
RFCTOR= (DJRN* DJRNH/DJR/DJR2N) *RAHSQ 
SGR2R=RFCTOR*GR2R (MXGR2M) 
SDGR2R=RFCTOR*GDR2R (MXGR2M) 
DO 2500 JR=MXGR2M-1,N+1,-1 
DJR=DJR+ONEN 
DJRN=DJRN+ONEN 
DJRNH=DJRNH+ONEN 
DJIR2ZN=DJR2N+ONEN 
RFCTOR= (DJRN* DJRNH/DJR/DJIR2N) *RAHSQ 
SGR2R= (GR2R(JR) -SGR2R) *RFCTOR 
SDGR2R= (GDR2R (JR) ~-SDGR2R) *RFCTOR 
2500 CONTINUE 
ho(N =. 50. 0) THEN 
GROUT= (GR2R (0) -SGR2R) *GRSF 
GDROUT= (GDR2R (0) -SDGR2R) *GRSF 
ELSE 
RFCTOR=HALF* (HALF-DN) *RAHSQ 
SGR2R= (GR2R(N) -SGR2R) *RFCTOR 
SDGR2R= (GDR2R(N) -SDGR2R) *RFCTOR 
DJR=DN 
DJRH=DJR~HALF 
DJRN=DJR+DN 
DJRNN=DJR~-DN 
BO 2600 JR=N-1,1,-1 
DJR=DJR+ONEN 
DJRH=DJRH+ONEN 
DJIRN=DJRN+ONEN 
DJIRNN=DJRNN+ONEN 
RFCTOR= (DJR*DJRH/DJRN/DJRNN) * RAHSQ 
SGR2R= (GR2R (JR) -SGR2R) *RFCTOR 
SDGR2R= (GDR2R (JR) -SDGR2R) *RFCTOR 
2600 CONTINUE 
GROUT= (GR2R (0) -SGR2R) *GRSF/DN 
GDROUT= (GDR2R (0) -SDGR2R) *GRSF/DN 
END IF 
END IF 
RETURN 
END 


Ik SUBROUTINE BCKSCFL 


SUBROUTINE BCKSCFL(IR,CESTH,CESPH) 
COR RRR RRR RRO III III UII III IORI III III OR II Ik 
(e 
INCLUDE ‘REALTP.INC' 
ENGLUDE “CMPXTP.INC* 
INCLUDE 'MAINDM.INC' 


REAL*4 AKPHPR,AKPHPI,AKZPR,AKZPI,ALPHPR, ALPHPI,ALZPR,ALZPI, 


oe 


100 


200 


300 
400 
Cc 


cS 


Cc 


+ 


AKPHNR, AKPHNI, AKZNR, AKZNI, ALPHNR, ALPHNI, ALZNR, ALZNI 


DIMENSION CIW(KCRNT) , CKLN (KCRNT) 
DIMENSION CIWO (KXCRT) , CKLNO (KXCRT) , CXPQNO (KXCRT,KXCRT) 
DIMENSION CRPQ1 (KCRNT,KCRNT) , CRPQ2 (KCRNT, KCRNT), 


+ 


CXPQN (KCRNT, KCRNT) , CKRPQ (KCRNT, KCRNT) 


DIMENSION -GRPHP (61,61) ,CKZP (61,61) (CLEHME (6liGie 


+ 
+ 


GCEZP (61.61), CKPHN (61,61) ;CKZN(61551)0 
CLPHN (61,61),CLZN(61, 61) 


COMMON /INPUT2/ IE,1IM, THETAI, THESINI, THECOSI,RTHEI 

COMMON /INPUT4/ IZ,1IK,IS,NYSM 

COMMON /NCONST/ DN, DNH,N 

COMMON /XPQTMP/ CXPON, CXRPQ, CXPOQNO 

COMMON /XPOTMP1/ CRPOQ1,CRPQ2 

COMMON /CKLMTX/ CKPHP, CKZP,CLPHP,CLZP, CKPHN, CKZN, CLPHN, CLZN 
COMMON /CRNTDM/ NMAX,MXNG, IQMAX, IQMAX1, IQMAX2, IXCRNT, ICRNT,MXQEG, 
+ MXQOG 

COMMON /INPUT3/ RTHE,RDELT,RPHI,RDELP,NTHTAO,NPHI, THESIN, THECOS, 
+ RHPI 

CESTH=CZERO 

CESPH=CZERO 

DOeLOO JO=1, KXERT 

CKLNO (JQ) =CZERO 

CONTINUE 

DO 200 JO=1, RERNT 

CKLN (JQ) =CZERO 

CONTINUE 

DO 400 IF=-30,30 


BO S00sd2——-307, 30 
CKPHP (IF, JZ) =CZERO 
CKZP (IF, JZ) =CZERO 
GEPHE (fe, JZ) =CZERO 
CLZP (IF,JZ)=CZERO 
CKPHN (IF, JZ) =CZERO 
CKZN(IF,JZ)=CZERO 
CLPHN (IF, JZ) =CZERO 
CLZN (IF,JZ) =CZERO 

CONTINUE 
CONTINUE 


iP (RTHEL EO. 0.D0) GO TO 2000 
C If incident angle is not equal to 0, use this loop 
DO 1800 NA=0,NMAX 


CALL INCIDNT (NA,CIWO,CIW) 


TP tna EO. 0) THEN 
CALL XPQ0 


ELSE 


CALL XPQN (NA) 

ENDIF 
C If the cylinder is made of perfect conductor and no coating on it 
TeeeGl Ze. BO, (0) THEN 


DN=NA 


C Use IMSL library to solve linear system 
CALL DLSACG(IXCRNT, CXPOQNO,KXCRT,CIWO,1,CKLNO) 
CALL ESCFAR(NA,CKLNO, CKLN, CETHN, CEPHN) 


Mo Bes 


If we calculate in X and Y components as theta approaches 0 or pi, 
then theta component is cosin phi in X direction plus sine phi in Y 


direction, 


cosin phi 


and phi component is negative Sine phi in X direction plus 
an Yodirection. 


Sree RTHE 2 2O.. 0..D0) “THEN 
CEGEN=COs(RPHI) “CETHN+SIN(CRPHI) *“CEPHN 


CEPHN=- 


SIN (RPHI) * CETHN+COS (RPHI ) *CEPHN 


ELGEEr (RI ne «nO. PI) THEN 
CETHN=-COS (RPHI) *CETHN-SIN(RPHI ) *CEPHN 
CEPHN=-SIN(RPHI) *CETHN+COS (RPHI) * CEPHN 

END IF 
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CEsTH=CESTH+CETHN 
CESPH=CESPH+CEPHN 
Le NAN. 0} THEN 
toc - 50. 1) THEN 
DG5600 L=2;,KXCRT, 2 
CKLNO (I) =-CKLNO (T) 
600 CONTINUE 
END IF 
LF AIM .EO.. 1) ° THEM 
DO 800 T= KXCRT, 2 
CKLNO (I) =-CKLNO (TI) 
800 CONTINUE 
BN@ iF 
NA1=-NA 
DN=NA1 
CALL ESCFAR(NA1, CKLNO, CKLN, CETHN, CEPHN) 
TP (RTHES. 20. 0.D0) THEN 
CETHN=COs (REEL) -“CEtaN+ sin (RPHI) *CEPHN 
CEPHN=-sthii(kent) “CEiMN+COs (RPHL) *CEPHN 
BuShLe s(RTHE .BO. Pi) THEN 
CETHN=—-COS (RPHI) *CETHN-SIN(RPHI) *CEPHN 
CEPHN=-SIN(RPHI) *CETHN+COS (RPHI) * CEPHN 
END iF 


CESTH=CESTH+CETHN 
CES PH=CESPH+CEPHN 
END Ge 
END IF 
Cc If the cylinder is coated with anisotropic material 
teetlZ eo. |} THEN 
DO 1100 I=1,KCRNT 
DO 1000 J=1,KCRNT 
CZRPO(l, J) =CXPON(1,J) +CRPOL (I,J) 


1000 CONTINUE 
1100 CONTINUE 
DN=NA 


CALL DLSACG (ICRNT, CXRPQ,KCRNT,CIW,1,CKLN) 
fee Dh ee bo, 1) .AND.. (IR .EQO: 45) )° THEN 
CALL KLCRNT (DN, CKLN,CIW) 
END IF 
CALL ESCFAR(NA,CKLNO, CKLN, CETHN, CEPHN) 
if (RTHE .EO. 0.DO) THEN 
GETHN=COS (RPHI) *CETHN+SIN(RPHI) *CEPHN 
GBPHN=—-SIN(RPH1) *“CETHN+COS (RPHI) *CEPHN 
Buosie (RTHE .£O.) Pi) Ss THEN 
CETHN=-CCS(RPHI) *CETHN-SIN(RPHI) *CEPHN 
CEPHN=-SIN(RPHI) * CETHN+COS ( RPHI) *CEPHN 
END IF 


CESTH=CESTH+CETHN 
CESPH=CESPH+CEPHN 
rr (NA -NE. 0) THEN 
IF (NYSM .EQ. 0) THEN 
beet 500 T=) KERN. 
DO 1200 J=1,KCRNT 
CXRPOCI, J) =CXPON (1 (J) +CRPO2Z(1,J) 
1200 CONTINUE 
1300 CONTINUE 
ELSE 
DO 72500 -1=1;7,KCRNT 
DO 1400 J=1,KCRNT 
CAREOC 1 7.) =CXPON(L, J )+CRPOL (I,J) 
1400 CONTINUE 
1500 CONTINUE 
END IF 
CALL DLSACG(ICRNT, CXRPQ, KCRNT,CIW,1,CKLN) 
IF (IE -.EQO. 1) 3@HEN 
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DOw~21600 1=2, KCRNIT 4 
Il=I+1 
CKLN (I) =-CKLN (TI) 
CKLN (I1) =-CKLN(I1) 
1600 CONTINUE 
END IF 
IF (IM .EQ. 1) THEN 
DO 1700 I=1,KCRNT,4 
PZ — is 
CKLN (I) =-CKLN (I) 
CKLN (I2 ) =-CKLN (I2) 
1700 CONTINUE 
END IF 
NA1=-NA 
DN=NA1 
Te eC 2 EO leer Darts. o> 2199).) THEN 
CALL KLCRNT (DN, CKLN,CIW) 
END IF 
CALL ESCFAR(NA1, CKLNO,CKLN, CETHN, CEPHN) 
Treetkine 60. 0.D0) THEN 
CETHN=COS (RPHI) *CETHN+SIN(RPHI) * CEPHN 
CEPHN=-SIN(RPHI) *CETHN+COS (RPHI) * CEPHN 
ELSEIF (RTHE .EQ. PI) THEN 
CETHN=-COS (RPHI) *CETHN~SIN (RPHI) *CEPHN 
CEPHN=-SIN(RPHI) *CETHN+COS (RPHI) *CEPHN 
END IF 


CESTH=CESTH+CETHN 
CES PH=CES PH+CEPHN 
END IF 
END IF 
1800 CONTINUE 
C If the cylinder is coated with anisotropic material and we want to calculate 
C equivalent currents K nad L on inner and outer surfaces. 
ci Oh. BO) a CMe 
IF (IS .EQ. 199) THEN 
OPEN (31, FILE='khzz41.dz', STATUS=' UNKNOWN ' ) 
OPEN (32, FILE='ihzz4l1.dz’ , STATUS=' UNKNOWN ' ) 
END IF 
DO 1900 Tr=-30,30 
DO 1850 JZ=-30,30 
AKPHPR=REAL (CKPHP(IF,JZ)) 
AKPHPI=IMAG (CKPHP(IF,JZ) ) 
AKZPR=REAL (CKZP(IF,JZ) ) 
AKZPI=IMAG (CKZP(IF,JZ) ) 
AKPHNR=REAL (CKPHN(IF,JdJZ) ) 
AKPHNI=IMAG (CKPHN(IF,JZ) ) 
AKZNR=REAL (CKZN(IF,JZ) ) 
AKZNI=IMAG (CKZN (IF, JZ) ) 
ALPHPR=REAL (CLPHP(IF,JZ) ) 
ALPHPI=IMAG (CLPHP(IF,JZ) ) 
ALZPR=REAL (CLZP(IF,JZ) ) 
ALZPI=IMAG (CLZP(IF,JZ) ) 
ALPHNR=REAL (CLPHN (IF,JZ) ) 
ALPHNI=IMAG (CLPHN(IF,JZ) ) 
ALZNR=REAL (CLZN(IF,JZ) ) 
ALZNI=IMAG (CLZN(IF,JZ) ) 


WRITE(31,*) AKPHPR, ' ', AKPHPL 
WRITES 1 *)) AKZPR, * ", AKZPI 
WRETE(S2, *) ALPHPR,* 7 AbPHEL 
WRETECS2),*) ALZPR, ' PA) sat 
WRITE (32, *) AKPHNE, ° ", AKPHNI 
WRITE(31,%*) AKZNR, ’ ', ,AKZNI 
WRITE(32,7) ALPHNR, ' ’ , ALPHNI 
WRITE<S 2") ALZNR, * ', ALZNI 
1650 CONTINUE 
1900 CONTINUE 
CLOSE (31) 


a5 


CLOSE (32) 
END IF 

RETURN 
© ££ incident angle is equal to zero, the program should go to this loop 
C because it need to compute when n=+1l and n=-1 only. 
2000 NA=1 

DN=NA 

CALL INCIDNT (NA, CIWO,CIW) 

CALL XPOQN (NA) 


LP (la. bO.. 0)) THEN 
CALL, DLSACG (IXCRNT, CXPONO, KXCRT, CIWw0, 1, CKLNO) 
CALL ESCFAR(NA, CKLNO, CKLN, CETHN, CEPHN) 


If we calculate in X and Y components as theta approaches 0 or pi, 
then theta component is cosin phi in X direction plus sine phi in Y 
@airection, and phi component is negative sine phi in X direction plus 
CcOSiN Phi in Y direction. 
IF (RTHE .EQ. 0.DO) THEN 
CETHN=COS (RPHI) *CETHN+SIN(RPHI) *CEPHN 
CEPHN=-SIN (RPHI) *CETHN+COS (RPHI) *CEPHN 
ELSEIF (RTHE .EQ. PI) THEN 
CETHN=-COS (RPHI) *CETHN-SIN(RPHI) *CEPHN 
CEPHN=-SIN(RPHI) *CETHN+COS (RPHI) *CEPHN 
END IF 


Hage Mat 


CESTH=CESTH+CETHN 
CESPH=CESPH+CEPHN 
IF (NA .NE. 0) THEN 
IF (IE .EQ. 1) THEN 
DOe2s00et=2- KXCRT 
CKLNO (I) =-CKLNO (TI) 
2400 CONTINUE 
END IF 
IF (IM .EQ. 1) THEN 
Dome OO sl=),KXCRT,2 
CKLNO (I) =-CKLNO (I) 
2500 CONTINUE 


END IF 
NAL=-NA 
DN=NA1 


CALL ESCFAR(NA1,CKLNO, CKLN, CETHN, CEPHN) 
LF SURTHE -5O 90550) o CHEN 

CETHN=COS (RPHE) *CETHN+SIN( RPHE) *CEPHN 

CEPHN=-SIN (RPHI) *CETHN+COS(RPHI) *CEPHN 
BLSELFE (RTHES £O.. FL) HEN 

CETHN=-COS(RPHI1) *CETHN-SIN(RPHI) *CEPHN 

CEPHN=-SIN(RPHI) *CETHN+COS(RPHI)*CEPHN 
END GIF 


CESTH=CESTH+CETHN 
CES PH=C ESPH+CEPHN 
END IF 

END IF 


TPZ .207.2). THEN 
DO 3100 I=1,KCRNT 
DO 3000 J=1,KCRNT 
CXRPOQ(I,J)=CXPON(1I,J)+CRPQ1 (I,J) 
3000 CONTINUE 
3100 CONTINUE 
DN=NA 
CALL DLSACG(ICRNT,CXRPQ,KCRNT, CIW,1,CKLN) 


PEM@erhee EO: 1) .AND.. (IR .EO. 45) )  0HEN 


CALL KLCRNT (DN, CKLN,CIW) 
END IF 
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3200 
3300 


3400 
5500 


3600 


3700 


CALL ESCFAR (NA, CKLNO, CKLN, CETHN, CEPHN) 
iP (RTHE 2. EO. 0-20), THen 

CETHN=COS (RPHI) * CETHN+SIN(RPHI) *CEPHN 

CEPHN=-SIN (RPHI) * CETHN+COS (RPHI) * CEPHN 
BPLSEIF (RTHE 2EO- el) THEN 

CETHN=-COS (RPHI) * CETHN-SIN (RPHI) *CEPHN 

CEPHN=-SIN (RPHI)* CETHN+COs (REH®) *CEPHN 
END oir 


CESTH=CESTH+CETHN 
CESPH=CESPH+CEPHN 
IF (NA .NE. 0) THEN 
IF (NYSM .EQ. 0) THEN 
DO eso0e 0-1 KCRNT 
DO 3200 J=1,KCRNT 
CXRPQ (I,J) =CXPON (I,J) +CRPQ2 (I,J) 
CONTINUE 
CONTINUE 
ELSE 
DO 3500 I=1,KCRNT 
DO 3400 J=1,KCRNT 
CXRPQ (I,J) =CXPON (I,J) +CRPQ1 (I,J) 
CONTINUE 
CONTINUE 
END IF 
CALL DLSACG(ICRNT,CXRPQ,KCRNT, CIW,1,CKLN) 
Tr ere). EO. 1) THEN 
DO 3600 I=2,KCRNT,4 
I1l=I+1 
CKLN (I) =-CKLN(T) 
CKLN(I1) =-CKLN(I1) 
CONTINUE 
END IF 
Le 40M .50:.. 1.) THEN 
DO. 3700 I=1,KCRNT,4 
EZ=1+2 
CKLN (I) =-CKLN (TI) 
CKLN (I2) =-CKLN(I2) 
CONTINUE 
END IF 
NA1=-NA 
DN=NA1 


Pre he el) AND oo (1S EO 199) THEN 
CALL KLCRNT (DN, CKLN, CIW) 
END IF 


CALL ESCFAR(NA1, CKLNO, CKLN, CETHN, CEPHN) 
IF (RTHE .EQ. 0.D0) THEN 

CETHN=COS (RPHI) * CETHN+SIN (RPHI) *CEPHN 

CEPoN=-Sin (hen) Ce lTHN+COs (RPH1) *CEPHN 
PECEIomcn tHe .60, PI) THEN 

GE THi=—COs (RPH1I)*CETHN-SIN(RPHI ) *CEPHN 

€EPHH——siN (RPHI) *CETHN+COS (RPHI) *CEPHN 
END IF 


CESTH=CESTH+CETHN 
CESPH=CESPH+CEPHN 
BND) IF 

END iF 


Pero. 1) THEN 
Pres 20. 199) THEN 
OPEN (31, FILE=‘khzz41.dz’ , STATUS='UNKNOWN' ) 
OPEN (32, FILE='1hzz41.dz', STATUS=' UNKNOWN' ) 
END IF 
DO 4200 IF=-30,30 
DO 4100 JZ=-30,30 


o/ 


AKPHPR=REAL (CKPHP(IF,JZ) ) 
AKPHPI=IMAG (CKPHP(IF,JZ) ) 
AKZPR=REAL (CKZP(IF,JZ) ) 
AKZPI=IMAG(CKZP(IF,JZ) ) 
AKPHNR=REAL (CKPHN(IF,JZ) ) 
AKPHNI=IMAG (CKPHN(IF,JZ) ) 
AKZNR=REAL (CKZN(IF,JZ) ) 
AKZNI=IMAG (CKZN(IF,JZ) ) 
ALPHPR=REAL (CLPHP(IF,JZ) ) 
ALPHPI=IMAG (CLPHP(IF,JZ) ) 
ALZPR=REAL(CLZP(IF,JZ) ) 
ALZPI=IMAG(ECLZP(IF, JZ) 
ALPHNR=REAL (CLPHN (IF, JZ) ) 
ALPHNI=IMAG (CLPHN(IF,JZ) ) 
ALZNR=REAL (CLZN (IF, JZ) ) 
ALZNI=IMAG (CLZN(IF,JZ) ) 


WRITE(31,*) AKPHER, * ‘7, AaReHEI 
WRITE (31%) AKZEPR, * 17 BR LEE 
WRITE (32,*) ALPHPR, ' ', ALPHPI 
WRETE(32; *)SALZPR; *  ALZEL 
WRITE(31,*) AKPHNR, ' ', AKPHNI 
WRITE(31,*) AKZNR, ' ', AKZNI 
WRITE (32;%~) ALPHNR 5° ', ALPHNI 
VRETEQS2, ~ ALZNR, *, ALZNI 
4100 CONTINUE 
4200 CONTINUE 
CLOSE (31) 
CLOSE(32) 
END? iF 
RETURN 
END 


Ce 
SUBROUTINE INCIDNT (NA, CIWO,CIW) 

CREEK K KKK RR RRR KK ERR KR KK RE KKK KR KEKE KK RK RRR KK KK KK KK RK KKK KR RK RR RK KR Kw 
C This subroutine sets up the incident wave on an object which is coated 
C with anisotropic material, or a perfect conductor. 
CE RR RRR KERKEKREKKEKKEEKKKEKKKKKKKEKKKEKKKKEKKKKKKKEKKKKEKEKEKKKKKKEKKEKKKKE 

INCLUDE ‘REALTP.INC' 

INCLUDE 'CMPXTP.INC' 

INCLUDE 'MAINDM.INC' 


Bs 
DIMENSION CIWO (KXCRT) ,CIW(KCRNT) , DJINS (KNDIM1+1) , DJPS (KODIM1+2) 
COMMON /INPUT2/ IE,IM, THETAI,THESINI, THECOSI,RTHEI 
COMMON /INPUT4/ IZ,IK,IS,NYSM 
COMMON /INPUT3/ RTHE,RDELT,RPHI,RDELP,NTHTAO, NPHI, THESIN, THECOS, 
+ RAPT 
COMMON /RTHETA/ DLICOSI,DL2SINI, DLL 
COMMON /GCONST/ DKH,DKA,HH, HA,HSQ, DASO,HSON, DASON, HESQ, HDASQ, 
~~ RAH , RAHSQ, DHA, DAH 
COMMON /CRNTDM/ NMAX,MXNG, IOMAX, IOQMAX1, IOMAX2, IXCRNT, ICRNT,MXOEG, 
ae MX0O0G 
ee 
NP=NA+1 
NP1=NP+1 
NM=NA-1 
S 
IF (IZ .EQ. 0) GO TO 4000 
& 
C This part computes the incident field on an anisotropic object 
C 
C Initialize the column matrix of sum current on the anisotropic 
DOLPOO. fJ=1, KCRNT 
CIW (IJ) =CZERO 
100 CONTINUE 
G 


C If the incident angle is 90 degree or 0 degree 
IF (RTHEI .EQ. ZERO) GO TO 3000 
CALL DBSJNS(DL2SINI, KNDIM1, DJNS) 
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DJN=DJINS (NP) 
IF (NA. EQ. O) THEN 
DDJN=-DJNS (NP1) 
ELSE 
DDJN=HALF* (DJNS (NA) -DJNS (NP1) ) 
ENDIF 
C Use IMSL library to compute Bessel'‘s function series. 
CALL DBSJNS (DL1ICOSI, IOQMAX2,DJPS) 
DO 600 IP=0,KQODIM 
NIP=IP+1 
DNIP=NIP 
NIPL=NIP+1 
NIP2=NIP+2 
N11=NA+IP-1 
N12=NA+IP 
IE1=MOD (N11, 4) 
TEZ=MOD(N1i2Z, 4} 


CPia=Gies SLL 
CP2=C1z2-"1E2 
IPW1=4*IP+1 
IPW2=IPW1+1 
IPW3=IPwW2+1 
IPW4=IPW3+1 


DJP=DJPS (NIP) 
DJP1=DJPS (NIP1) 
DJP2=DJPS (NIP2) 
teeth .bO.el) THEN 
CEEPHI=CF1* (DJP+DJP2) *DDJN/DNIP 
CEHPHI=-TWO*NA*CF2*DJP1* DJN/DLL 
CEHZ=-CF2* THESINI* (DJP+DJP2) *DJN/DNIP 
ELSE 
CEEPHI=CZERO 
CEHPHI=CZERO 
CEHZ=CZERO 
END iF 


ip SiMe EO. 1) THEN 
CMEPHI=TWO*NA*CF2*DJP1*DJN/ DLL 
CMEZ=CF2* THESINI* (DJP+DJP2) *DJN/DNIP 
CMHPHI=CF1* (DJP+DJP2) *DDJN/DNIP 
ELSE 
CME PHI=CZERO 
CMEZ=CZERO 
CMHPHI=CZERO 
END IF 


CIW (IPW1) =TWO* (CEEPHI+ *CMEPHI ) 
CIW (IPW2) =TWO* CMEZ 
CIW (IPW3 ) =TWO* (CEHPHI+CMHPHIT ) 
CIW(IPW4) =TWO*CEHZ 

600 CONTINUE 
RETURN 


3000 CALL DBSJNS (DKH, IQMAX2,DJPS) 
DO 3200 IP=0,KQDIM 
JP=IP+2 
DJP=DJPS (JP) 
IE11=MOD (IP, 4) 
IE12=MOD(IP+1, 4) 
Chil=-Cil2**lB11 
GFrZ2=C12**1E1Z 
IPW1=4*IP+1 
IPW3=IPW1+2 


TP (1B .EO. 2) THEN 
CEEPHI=CF1*DJP/DKH 
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CEHPHI=-CF2*DJP/DKH 
ELSE 
CEEPHI=CZERO 
CEHPHI=CZERO 
END IF 


ie CiM .FOs 1) THEN 
CMEPHI=CF2* DJP/DKH 
CMHEHI=Crl “Dar DAH 

ELSE 
CMEPHI=CZERO 
CMHPHI=CZERO 

END iF 


CIW (IPW1) =TWO* (CEEPHI+CMEPHI) 
CIW (IPW3 ) =TWO* (CEHPHI+CMHPHTI) 
3200 CONTINUE 
RETURN 
ie 
C This part compute incident fields on a perfect conductor 
C Initialize the column matrix of sum current on the conductor 
4000 9 D684200 ix = 1 RxACRT 
CIWO (IX) =CZERO 
4200 CONTINUE 
& 
Cc If the incident angle is 90 degree or 0 degree 
IF (RTHEI .EQ. ZERO) GO TO 6000 
C 
CALL DBSJNS(DL2SINI, MXNG+1, DJNS) 
DIN=DJINS (NP) 
IF (NA .EQ. 0) THEN 
DDJN=-DUNS (NP1) 
ELSE 
DDJN=HALF* (DJNS (NA) -DJNS(NP1) ) 
ENDIF 


GALE, DESINS (DEIECOSI, IOMAX2, DIPS) 
DO 4400 IP=0,KQDIM 

NEPaIPs |) 

DNIP=NIP 

NIP1=NIP+1 

NIP2Z=NIP+2 

N11=NA+IP-1 

N12=NA+IP 

IE1=MOD (N11, 4) 

IE2=MOD(N12, 4) 


Cll=Ci 2 “ere 
CF2=CI2**fE2 
TPW1=2* IP+1 
IPW2=IPW1+1 


DJP=DJIPS (NIP) 
DJP1=DIPS (NIP1) 
DIP2Z2=DIPS (NIP2) 
TE (TE 2EO. sl) THEN 
CEEPHI=CF1* (DJP+DJP2) *DDJN/DNIP 
ELSE 
CEEPHI=CZERO 
END IF 


ie CEM. EO. 1) THEN 
CMEPHI=TWO*NA* CF2*DJP1*DJN/DLL 
CMEZ=CF2*THESINI* (DJP+DJP2) *DIN/DNIP 
ELSE 
CMEPHI=CZERO 
CMEZ=CZERO 
END IF 


100 


CIWO (I PW1) =TWO* (CEEPHI+CMEPHI) 
CIWO (I PW2 ) =TWO*CMEZ 

4400 CONTINUE 
RETURN 


6000 CALL DBSJNS (DKH, IQMAX2,DJPS) 
DO 6200 IP=0,KQDIM 
JP=IP+Z 
DJP=DJPS (JP) 
IE11=MOD (IP, 4) 
IE12=MOD (IP+1, 4) 
CPl=6Gi2> Ter 
CF2=Ci2*> 71212 
IPW1=2*IP+1 


& 
PP eto es fo. 2) THEN 
CEEPHI=CF1*DJP/DKH 
ELSE 
CEEPHI=CZERO 
END IF 
< 
Cecio... 1) THEN 
CMEPHI=CF2*DJP/DKH 
ELSE 
CMEPHI=CZERO 
END LF 
c 


CIWO (I PW1) =TWO* (CEEPHI+CMEPHI ) 
6200 CONTINUE 
RETURN 
END 
= 
SUBROUTINE XPQO0 
CER EKRKKE KK KK KK EKER KK KEKE KKK KEKE KKK KKK KK KK KKK KKK KKK KKK KEK KKK KKK KK kK 
C This subroutine computes the matrix XN(P,Q) for N = 0 following a 
C call to XPQINI. This matrix is kept in the common block: 
Ec COMMON /XPOQOTMP/ CXPON 
CERKKKKKKKK KKK KK KEK KK KKK KKK KKK KKK KK KKK KEK KKK KKK KK KKK KKK KKK KKK KKK KKK KKK KK KK 
INCLUDE ‘REALTP.INC' 
INCLUDE ‘CMPXTP.INC'’ 
INCLUDE 'MAINDM.INC' 


© 
DIMENSION CXPON(KCRNT,KCRNT) , CXRPQ (KCRNT, KCRNT) 
DIMENSION CXPQNO (KXCRT,KXCRT) 
DIMENSION CGNE(0:MAXPEG,0:MAXPEG,KNDIM1+1), 
+ CGNO(0:MAXPOG, 0: MAXPOG, KNDIM1+1) 
DIMENSION CDGNE(0:MAXPEG,0:MAXPEG, KNDIM1+1), 
af CDGNO(0:MAXPOG, 0:MAXPOG, KNDIM1+1) 
COMMON /GPQOTMP/ CGNE, CDGNE,CGNO,CDGNO 
COMMON /XPQOTMP/ CXPON, CXRPQ, CXPQNO 
COMMON /CRNTDM/ NMAX,MXNG, IQMAX, IOQMAX1, IOQMAX2, IXCRNT, ICRNT, MXQEG, 
ae MXQO0G 
COMMON /GCONST/ DKH, DKA, HH,HA,HSQ, DASQ,HSQON, DASON, HHSQ, HDASQ, 
+ RAH, RAHSQ, DHA, DAH 
COMMON /INPUT4/ IZ,1IK,1IS,NYSM 
SAVE /XPQTMP/,/GPQTMP/,/GCONST/ 
& 
IF (N .NE. 0) THEN 
WRITH(*) Je Input N is not equal to 0 ian XPO0O- 
WRITE(*,*) ‘Execution is stopped. ' 
STOP 
END IF 
c 
IF “(iZ -5EO. 0) Go To 4000 
c 
C Initialize the XN(P,Q) matrix. 
C Null terms will be skipped later. 
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DO 200 IQ=1,KCRNT 
por Loo 1P=1,KCRNT 
GxXPOn tir, £0) =CZERO 


100 CONTINUE 
200 CONTINUE 
C 


C Form the XN(P,Q) matrix for n = 0 and using on a perfect conductor 
CF31=DKA*CI1 
ce 
DO 1300 IQE=0,MXQEG-1 
IQE1=IQE+1 
LTO=2* 108 
IQX1=4*IO+1 
IQOX2=IOX1+1 
IQX3=IQOX2+1 
IQOX4=I0X3+1 
DOL=10+1 
FQ22=D01/HHSO 
DO 1100 IPE=0,MXQEG-1 
IPE1=IPE+1 
LP=2* 1 PE 
IPX1=4*IP+1 
IPX2=IPX1+1 
IPX3=IPX2+1 
IPX4=IPX3+1 
DP1=IP+1 
F41=DKH/DP1 
CF11=F41*CF31 
FS1=QUAR*F41*DKA 
CXPON (IPX1, IQX1)=(CGNE (IPE1,IQE,2) -CGNE(IPE,IQE,2))*CF11 
CXPOQN (IPX4, IQX1)=(CGNE(IPE1,IQE,2)-CGNE(IPE, IQE,2)+HA* ( 


+ CDGNE (IPE1, IQE, 2) -CDGNE (IPE, IQE,2)))*F41 
CXPON (1PX2 , 10K2) = (TWO"EO227DP1*CGNO (IPE, IQE, 1) +HALF* ( 

+ CGNE (IPE, IQE1,1)+CGNE(IPE1,IQE,1)- 

+ CGNE (T PE], IQOEL, U)-CGNE (1 PE, DOE, 2) )) ~cr lt 
CXPON (I PX3, 10X2Z)=(CDGNE (IPE, IOE,1)+CDGNE( IPE] aOE1, 1)=— 

+ CDGNE(TPEI, LOE, 1)-CDGNE (IPE, IOET 71) ) bo 


CXPON (I PX2 , 1QX3) =-CXPQN (IPX4, 1QX1) 
CXPOQN (IPX3 , 1QX3) =CKPQN (IPX1, IQX1) 
CXPOQN (IPX1, IQX4) =-CXPQN (IPX3, IQX2) 
CXPOQN (IPX4, IQX4) =CXPOQN (IPX2, IQX2) 
1100 CONTINUE 
1300 CONTINUE 
DO 2300 IQ0=0,MXQ0G-1 
IQO01=IQ0+1 
IQO=2*IO00+1 
IQX1=4*I0+1 
TOX2=I1OX1L+1 
IQX3=IQX2+1 
IQX4=I1QOX3+1 
DQO1=IQ+1 
FQ22=DQ01/HHSQ 
DO 2200 IPO=0,MXQ0G-1 
LEOlL=LFOr EL 
IP=2* IPO+1 
Eexi 42 1 P+i 
IPX2=IPX1+1 
IPX3=IPX2+1 
IPX4=IPX3+1 
DP1l=IP+1 
F41=DKH/DP1 
CF11=F41*CF31 
F51=QUAR*F41*DKA 
CXPON (1PX1, 1OX1) =(CGNO(1PO1, 100, 2)-CGNO( IPO, 100, 2) )2CrrE 
CXPON (IPX4,1I0QX1) = (CGNO (IPO1,IQ0, 2) -CGNO (IPO, IQ0, 2) +HA* ( 
= CDGNO(1.PO1, 100, 2) -CDGNO (LPO, 1907 2))) e241 
CXPON (IPX2, IQX2) = (TWO* FQ22*DP1*CGNE (IPO1, IQ01,1)+HALF™ ¢ 


102 


+ CGNO(IPO, 1001, 1) +CGNO({TPOl 1Ge7 1) 


+ CGNO(IPO1, 1001.,1) -CGNO(IPO71GG 7 a eri 
CXPON (IPX3 , I1QOX2) = (CDGNO (IPO, 100, 1) +CDGNO(IPO1,1I9Q01,1)- 
+ CDGNO(IPO1, 100, 1) -CDGNO( 1 PO, foe ee ow 


CXPOQN (IPX2, IOX3) =-CXPON (IPX4,10QX1) 
CXPON (IPX3 , IOX3 ) =CXPON (IPX1,IOX1) 
CXPON (IPX1, 10X4) =-CXPOQN (IPX3 , 10X2) 
CXPON (IPX4, IOX4) =CXPON (IPX2,IOX2) 


2200 CONTINUE 
2300 CONTINUE 
RETURN 


C Initialize the XN(P,Q) matrix. 
C Null terms will be skipped later. 
C 
4000 DO 4200 I0=1,KXCRT 
DO 4100 IP=1,KXCRT 
CXPONO (IP, 10) =CZERO 


4100 CONTINUE 
4200 CONTINUE 
ee 


C Form the XN(P,Q) matrix for n = 0 and using on an anisotropic coat. 
Ceol —pKA*erl 
e 
DO 4500 IQE=0,MXQEG-1 
IQE1=IQE+1 
TO=Z2“1OE 
TOR V=2* LO 
IOX2=I0OX1+1 
DO1=I0+1 
FQ22=D01/HHSO 
DO 4400 IPE=0,MXQEG-1 
IPE1=IPE+1 
IP=2* IPE 
IPX1=2* IP+1 
IPX2=IPX1+1 
DP1=IP+1 
F41=DKH/DP1 
CF11=F41*CF31 
F51=QUAR*F41*DKA 
CXPONO (IPX1, IQX1) =(CGNE(IPE1,IQE, 2) -CGNE(IPE, 1IQE,2))*CF11 
CXPONO (IPX2, IQX2) = (TWO* FO22*DP1*CGNO (IPE, IQE, 1) +HALF* ( 
+ CGNE (IPE, IQE1,1)+CGNE(IPE1,IQE,1)- 
+ CGNE(IPE1, IQE1,1)-CGNE(IPE,IQE,1)))*CF11 
4400 CONTINUE 
4500 CONTINUE 
DO 5300 IQ0=0,MXQ0G-1 
I901=1I00+1 
EQ=2* 100+1 
IQOX1=2*I0+1 
IOX2=I10X1+1 
DO1=I90+1 
FQ22=D01/HHSOQ 
DO 5200 IPO=0,MX0Q0G-1 
IPO1=IPO+1 
IP=2* IPO+1 
IPX1=2* IP+1 
IPX2=IPX1+1 
DP1=IP+1 
F41=DKH/DP1 
CF11=F41*CF31 
F51=QUAR* F41*DKA 
CXPONO (IT PX1, I0X1) =(CGNO(IPO1,1IQO, 2) -CGNO(1IPO, 100, 2) ) *CF1l 
CXPONO (IPX2, IQOX2) = (TWO* FQ22*DP1*CGNE (IPO1,10Q01,1)+HALF* ( 


+ CGNO (IPO, IQ01,1)+CGNO(IPO1,1Q0,1)- 

7 CGNO(IPO1, 1001, 1) -CGNO(1TPO, 10071)9)) eral 
5200 CONTINUE 
9300 CONTINUE 

RETURN 
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END 
es 

SUBROUTINE XPON(NIN) 
CRE KKK KKK KK KK Ke eR eR mk kk ee ee eK ee ke eR kk ee ek ek ke ek ee ek ke ee 
C This subroutine calls the subroutine GDN to update G(P,Q,N) for 
CyyN = NIN+LD, then forms the matrix XN(P,O) for N = NIN = 02 This matrix 
C is kept in the common block: 
ce COMMON /XPQOTMP/ CXPON 


CER KK KKK KEKE K KKK KKK KKK KKEK KKK KKK KK KA K KK KKK KK KKK KKK KK KKK KKK KKK KKK 


C WARNING: 

C IT IS ASSUMED THAT XPQINI AND XPQN FOR N FROM 1 TO NIN-~-1 HAVE BEEN 
C CALLED SO THAT G(P,0,N) AND ITS DERIVATIVE FOR N=NIN-1 AND N=NIN 

C HAVE BEEN STORED IN THE COMMON BLOCK /GPQTMP/ WITH PROPER N-INDICES. 


TNGCEU DE ~REALTP. INC ' 
TNEEUDE “CMPATP.INC* 
INCLUDE '‘MAINDM.INC' 


DIMENSION CXPON (KCRNT,KCRNT) , CKRPO (KCRNT, KCRNT) 

DIMENSION CXPONO (KXCRT,KXCRT) 

DIMENSION CGNE(0:MAXPEG,0:MAXPEG, KNDIM1+1), 

+ CGNO (0: MAXPOG, 0:MAXPOG, KNDIM1+1) 

DIMENSION CDGNE(0:MAXPEG,0:MAXPEG,KNDIM1+1), 

+ CDGNO (0: MAXPOG,0:MAXPOG, KNDIM1+1) 

COMMON /GPQTMP/ CGNE, CDGNE,CGNO, CDGNO 

COMMON /XPQTMP/ CXPOQN, CXRPQ, CXPONO 

COMMON /CRNTDM/ NMAX, MXNG, IOMAX, IOMAX1, IOMAX2, IXCRNT, ICRNT,MXOEG, 
+ MX0OO0OG 

COMMON /GCONST/ DKH, DKA,HH,HA,HSQ, DASO,HSON, DASON, HHSO, HDASO, 
+ RAH, RAHSO, DHA, DAH 

COMMON /NCONST/ DN, DNH,N 

COMMON /INPUT4/ IZ,IK,IS,NYSM 

SAVE /XPQOTMP/,/GPQTMP/,/GCONST/ 


N=NIN 

NA=ABS (N) 

Le {NA .LT. 1) THEN 

WRITE(*,*) ‘Input ABS(N) is less than 1 in XPOQN.' 
WRITE (*,*) ‘Execution is stopped.' 

STOP 

END IF 


NOT=NA+1 
NPI=NOI+1 
NMI=NA 


DNO=NA 
DNP=NOT 
DNM=NA-1 
TP eCGZ, EO. 0) GO TO 4000 


Initialize the XN(P,Q) matrix. 
Null terms will be skipped later. 


(C1 @ ll @ ag 


DO 800 IO=1,KCRNT 
Dow 7O0 LP=],KCRNT 

CXPON (IP, IQ) =CZERO 

ZOO CONTINUE 

800 CONTINUE 

c 

C Form the XN(P,Q) matrix. 
CES I=DKA*CIil 
F21=TWO*DNO 
FN11=F21*DNO/DASOQ 


DO 1300 IQE=0,MXQEG-1 
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1100 


1200 
1300 


ae 
+} 


+ 


=f. 


+ 


a 
ft 


+ 


IQOET=T1OB+1 
IQ=2* IQ0E 
IQX1=4*I0+1 
TOX2=I90X1+1 
IOX3=IOX2+1 
IOX4=IOX3+1 
DQO1=I90+1 
FO12=DN0 *DO1 
FQ22=D01/HHSO 
DO 1100 IPE=0,MXQEG-1 
IPE1=IPE+1 
IP=2* IPE 
IPX1=4*IP+1 
IPX2=IPX1+1 
IPX3=IPX2+1 
IPX4=IPX3+1 
DP1=IP+1 
F41=HH/DP1 
CFII=F41*CF31 
F51=HALF* DKA*F41 
CXPON (IPX1, IQX1)=( (CGNE(IPE, IQE,NOI)-CGNE(IPE1,IQE,NOI) ) *FN11+ 
CGNE (IPE1, IQE, NMI) -CGNE (IPE, IQE, NMI) + 
CGNE (IPE1, IQE,NPI) -CGNE(IPE, IOQE,NPI))*CF11 
CXPON (IPX4, IQOX1)=( (CGNE(IPE, IQE, NMI) -CGNE(IPE1, IQE,NMI) ) *DNM+ 
(CGNE(IPE1, IQE,NPI) -CGNE (IPE, IQE,NPI) ) *DNP+ 
HA* (CDGNE (IPE1, IQE,NMI) -CDGNE (IPE, IQE, NMI) + 
CDGNE (IPE1,IQE,NPI)-CDGNE(IPE,IQE,NPI)))*F41 
CXPON (IPX2,IQX2)=(FQ22*DP1*CGNO(IPE,IQE,NOI)+ 
CGNE (IPE, IQE1,N0I)+CGNE(IPE1,IQE,NOI)- 
CGNE (IPE1,IQE1,NOI) -CGNE(IPE,IQE,NOI))*CF11 
CXPON (IPX3, 10X2) = (CDGNE (IPE, IQE,NO1I)+CDGNE(IPE1,IQE1,NOI)- 
CDGNE (IPE1,IQE,NOI) -CDGNE(IPE,IQE1,NOI))*F51 
CXPON (IPX2, 10X3) =-CXPON (IPX4, 1I0X1) 
CXPON (IPX3, IOX3 ) =CXPON (IPX1, IQX1) 
CXPON (IPX1,1IQX4)=-CXPON(IPX3,IQX2) 
CXPON (IPX4, I0X4) =CXPON (IPX2, IOX2) 
CONTINUE 
DO 1200 IPO=0,MXQ0G-1 
IPO1=IPO+1 
IP=2*IPO+1 
IPX1=4*IP+1 
IPX2=IPX1+1 
IPX3=IPX2+1 
IPX4=IPX3+1 
DP1=IP+1 
F12=FQ12/DP1 
CXPON (IPX2,10X1) =F21*CGNE(ZIPO1,IQE,NOI) 
CXPON (IPX3, 1IQXK1)=(CGNE(IPO1,IQE,NMI)-CGNE(IPO1,IQE,NPI))*CF31 
CXPON (IPX1,1IQX2) =(CGNO(IPO1,IQE,NOI)-CGNO(IPO,IQE,NOI))*F12 
CXPON (IPX1,IQX3)=-CXPOQN(IPX3,IQOX1) 
CXPON (IPX4, I10X3) =CXPON (IPX2,IQX1) 
CXPON (IPX3 , 10X4) =CXPON (IPX1, IQX2 ) 
CONTINUE 
CONTINUE 
DO 2300 I9Q0=0,MX00G-1 
IQ001=1I1900+1 
TO=2* FO0+1 
IOX1=4*I0+1 
IOX2=I0X1+1 
IQX3=I0OX2+1 
TOX4=I0X3+1 
DOQ1=I90+1 
FQ12=DNO*DQO1 
FQO22=D01/HHSOQ 
DO 2100 IPE=0,MX0OEG-1 
IPE1=IPE+1 
IP=2* IPE 
IPX1=4*IP+1 
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IPX2=IPX1+1 
IPX3=IPX2+1 
IPX4=IPX3+1 
DP1=IP+1 
F12=FQ12/DP1 
CXPON (IPX2, IQX1) =F21*CGNO(IPE,IQ0O,NOTI) 
CXPON (IPX3, IQX1) =(CGNO(IPE, I00, NMI) -CGNO(IPE,IQO,NPI))*CF31 
CXPON (IPX1, IQX2) =(CGNE (IPE1, IQ01,NO0OI)-CGNE(IPE,IQ01,NOI))*F12 
CXPON (IPX1,IQX3) =-CXPOQN (IPX3, IQOX1) 
CXPON (IPX4, IQX3) =CXPON (IPX2,IQX1) 
CXPON (IPX3 , IOX4) =CXPOQN (IPX1, IQX2) 
2100 CONTINUE 
DO 2200 IPO=0,MXQ0OG-1 
IPOL=IPO+1 
IP=2*IPO+1 
IPX1=4*IP+1 
EPRZ=tPxi+1 
IPX3=IPX2+1 
IPX4=IPX3+1 
DP1=IP+1 
F41=HH/DP1 
Ceilerai*ersl 
FS1=HALF* DKA*F41 
CXPON (IPX1, IQX1)=( (CGNO (IPO, IQ0, NOI) -CGNO(IPO1, IQ0,NOI) ) *FN11+ 


+ CGNO (IPO1, IQ0, NMI) -CGNO (IPO, IQ0, NMI) + 

~ CGNO (EPG oo, NPI) -CGNO(LEO; 100, NPL} =Cr it 
CXPOQN (IPX4, IQX1) =( (CGNO (IPO, IQ0, NMI) ~CGNO(IPO1,1IQO,NMTI) ) *DNM+ 

+ (CGNO(TPOL, TOO, NPI) -CGNO( TPO, 100, NPT) ) -PNP= 

+ HA* (CDGNO (IPO1, IQ0,NMI) -CDGNO (IPO, IQ0, NMI) + 

+ CDGNO (IPO1, IQ0,NPI) -CDGNO (IPO, I100,NPI)))*F41 
CXPON (I PX2, ITOX2) =(FO22* DEL *CGNE (TPOL, TOOL, NOL} + 

+ CGNO (IPO, IQ01,N0I)+CGNO(IPO1, IQ0,NO0T) - 

+ CGNO(IPO1,IQ001,N0I)~-~CGNO(IPO, IQ0,NO0I))*CF11 
CXPOQN (IPX3 , IQX2) = (CDGNO (IPO, IQ0, NOI) +CDGNO(IPO1,IQ01,NO0I)- 

> CDGNO(TPOl, 100, NOI) -CDGNO( EPO; LOOP NGT) )+FSt 


CXPOQN (IPX2, IQX3) =-CXPOQN (IPX4,IQX1) 
CXPOQN (IPX3, IQX3 ) =CXPOQN (IPX1, IQX1) 
CXPQN (IPX1, IQX4) =-CXPOQN (IPX3 , IQX2) 
CXPOQN (IPX4, IQX4) =CXPOQN (IPX2, IQX2) 


2200 CONTINUE 
2300 CONTINUE 
RETURN 


C Initialize the XN(P,Q) matrix. 
C Null terms will be skipped later. 
c 
4000 DO 4800 IQ=1,KXCRT 
DO 4700 IP=1,KXCRT 
CXPONO (IP, IQ) =CZERO 


4700 CONTINUE 
4800 CONTINUE 
C 


C Form the XN(P,Q) matrix. 
GES L=DKA*Cil 
F21=TWO*DNO 
FN11=F21*DNO/DASQO 


DO 5300 IQE=0,MXQEG-1 
IQE1=IQE+1 
TO=2*10E 
TOx1=2 *10O+1 
IQOX2=IQOX1+1 
DQ1=I0+1 
FQ12=DN0*DO1 
FQ22=DQ01/HHSQ 

DO 5100 IPE=0,MXQEG-1 
IPEL=IPE+1 
IP=2* IPE 
IPX1=2* IP+1 
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5100 


5200 
5300 


6100 


IPX2=IPX1+1 
DP1=IP+1 
F41=HH/DP1 
CPI =F4i7 crs 
FS1=HALF* DKA*F41 
CXPOQNO (IPX1,1IQX1)=((CGNE(IPE,IQE,NOI)-CGNE(IPE1, IQE,NOI)) *FN11+ 
CGNE(IPE1,IQE,NMI) -CGNE(IPE, IQE,NMI) + 
CGNE (IPE1,IQE,NPI) -CGNE(IPE,IQE,NPI) ) *CF11 
CXPQNO (IPX2, IQ0X2) = (FQ22*DP1*CGNO (IPE, IQE,NOI) + 
CGNE (IPE, IQE1,N0I)+CGNE(IPE1,IQE,NOI)- 
CGNE(IPE1,IQE1,N0I)-CGNE(IPE, IQE,NOI))*CF11 
CONTINUE 
DO 5200 IPO=0,MXQ0G-1 
IPO1=IPO+1 
IP=2* IPO+1 
IPX1=2*IP+1 
IPX2=IPX1+1 
DP1=IP+1 
F12=FQ12/DP1 
CXPONO (IPX2, IQK1) =F21*CGNE(IPO1,IQE,NOI) 
CXPQNO (IPX1,IQX2)=(CGNO(IPO1,IQE,NO0OI) -CGNO(IPO, IQE,NOI) )*F12 
CONTINUE 
CONTINUE 
DO 6300 IQ0=0,MXQ0G-1 
IQO01=IQ0+1 
TO=2*100+1 
TOX1=2*IO+1 
IQX2=IQK1+1 
DQ1=I0+1 
FQ12=DNO0*DQ1 
FQ22=D01/HHSQO 
DO 6100 IPE=0,MXQEG-1 
IPE1=IPE+1 
TP=2"2PE 
IPX1=2*IP+1 
IPX2=IPX1+1 
DP1=IP+1 
FI2Z=FOl12/ DP 
CXPQNO (IPX2, IQX1) =F21*CGNO(IPE,1IQ0O,NOI) 
CXPQNO (IPX1, IQX2) = (CGNE(IPE1,1IQ01,NO0I) -CGNE(IPE,1IQ01,NO0I))*F12 
CONTINUE 
DO 6200 IPO=0,MXQ0G-1 
IPO1=IPO+1 
IP=2* IPO+1 
IPX1=2*IP+1 
IPX2=IPX1+1 
DP1=IP+1 
F41=HH/DP1 
CF1I1=F41*CF31 
FS1=HALF* DKA*F41 
CXPONO (IPX1, 10X1) =( (CGNO(IPO,IQ0,NOI) -CGNO(IPO1,IQ0O,NOI1)) *FN11+ 
CGNO (IPO1,10Q0,NMI) -CGNO (IPO, IQ0, NMI) + 
CGNO (IPO1, IQ0,NPI)-CGNO (IPO, IQ0,NPI)) *CF11 
CXPONO (IPX2, IQ0X2) = (FQ22*DP1*CGNE (IPO1,1I001,NO0I)+ 
CGNO (IPO, IQ01,N0I)+CGNO(IPO1,IQ0,NOI)- 
CGNO (IPO1,1I9001,NO0I) -CGNO (IPO, IQ0,NOI)) *CF11 
CONTINUE 
CONTINUE 
RETURN 
END 


SUBROUTINE ESCFAR (NIN, CKLNO, CKLN, CETHN, CEPHN) 


CEE RR ERK KKK KKK KK KEKE KEKE KKK KEKE KEKE RE KEKE KEKE KEKE KKK KEKE KKK ERE KEKE KKK KEKE KEKE KKK KKK KKK KKK KK EK K 


C This subroutine computes the scattered fields in far zone 


@ 


INCLUDE ‘REALTP.INC' 
INCLUDE 'CMPXTP.INC' 
INCLUDE 'MAINDM.INC' 
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100 


DIMENSION CKLNO (KXCRT) 

DIMENSION CKLN(KCRNT) , DINS (KNDIM1) , DJPS (KQDIM1+2) 

COMMON /CRNTDM/ NMAX,MXNG, IQMAX, IOMAX1, IOQMAX2, IXCRNT, ICRNT,MXQEG, 
MXQOG 

COMMON /GCONST/ DKH,DKA,HH,HA,HSQ, DASQ, HSQN, DASQN, HHSQ, HDASQ, 

RAH, RAHSQ, DHA, DAH 

COMMON /INPUT3/ RTHE,RDELT,RPHI,RDELP, NTHTAO,NPHI, THESIN, THECOS, 
Ree 

COMMON /INPUT4/ 1I2Z,1K,IS,NYSM 


N=NIN 

DNi=N 

NA=ABS (NIN) 

NP=NA+1 

NP1=NP+1 

PHI1=RPHI*N 

Gi2Z2N=Ci2**N 

CENP=CONE*COS (PHI1)+C1I1*SIN(PHI1) 
CETHN=CZERO 

CEPHN=CZERO 


Pre GORTHE 2BO 0.) SOR (RE nE. bo.) Pil) )) GO. lomo ug 


TEOT=THECOS/ THESIN 

DLZSIN=DKA* THESIN 

DL1COS=DKH* THECOS 
DNL2=DN1*TCOT/ DKA 

CALL DBSJNS(DL2SIN, KNDIM1, DJNS) 


DIN=DJUNS (NP) 

IF (N .EQ. 0) THEN 
DDJN=-DJINS (NP1) 

ELSE 
DDJN=HALF* (DJNS (NA) -DJNS (NP1) ) 

ENDIF 


KN=NA/2 
KNP=NP/2 
TeevViNe bi. Oy ~ANDSR(KNE Se Gr. KN) ) 2THen 
DJN=-DJUN 
DDJN=-DDJN 
ENDIF 


CALL DBSJNS(BLICOS, KODIMI+2, DJPS) 


TE SCE? . EO eae 
DO 100 IP=0,IQMAX 
INP=IP+1 
NP2=INP+2 
N11=IP 
N12=IP+2 
IE1L=MOD (Nil, 4) 
CPi=Cr2** Tel) 
IPW1=2*IP+1 
IPW2=IPW1+1 
DJP=DJPS (INP) 
DJP2=DJPS (NP2) 
DJPH=HALF* (DJP+DJP2) 
CETHN=CETHN+CI2N* DHA* CENP*CI1*CF1* DJN* (DNL2 *CKLNO (IPW1) *DJP 

-THESIN*CKLNO (IPW2) *DJPH) 

CEPHN=CEPHN-CI2N* DHA* CENP*DDJN* CKLNO (IPW1) *CF1*DJP 
CONTINUE 

ELSE 
DO 200 IP=0, IOMAX 
INP=IP+1 
NP2=INP+2 
N1ii=IP 
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200 


2000 


2100 


NIZ=1P+2 

IE1=MOD(N11,4) 

CEi-—Ciz =~ ite 

TPw1l=4*IP+1 

IPW2=IPW1+1 

IPW3=IPW2+1 

IPW4=IPW3+1 

DJP=DJPS (INP) 

DJP2=DIPS (NP2) 

DJPH=HALF* (DJP+DJP2) 

CETHN=CETHN+CI2N* DHA* CENP* ( (CI1*DJN* DNL2* CKLN ( IPW1) -DDJN 
*CKLN (IPW3) ) *DJP-CI1*DJN* THESIN* CKLN (IPW2) *DJPH) *CF1 

CEPHN=CEPHN-CI2N*DHA* CENP* ( (DDJN*CKLN (IPW1) +CI1*DJN* DNL2 
*CKLN (IPW3) ) *~DJP-CI1* DJN* THESIN* CKLN (IPW4) *DJPH) *CF1 

CONTINUE 

END IF 
RETURN 


2F  (NASeNE:. 1) THEN 
GO TO 3000 
END IF 
CALL DBSJNS (DKH, KQDIM1, DJPS) 


LF (12 2 EO. 70) a0 eEN 
DO 2100 IP=0, IQMAX 
tINP=iP+1 
IPWl=2*IP+1 
JP=MOD(IP, 4) 
Chi=CLicrTJe 
CF2=CI2**JP 


DJP=DJIPS (INP) 
iE {RTRE 20-0.) THEN 
TE OCNeear. 0) THEN 
CETHN=CETHN-DAH* DJP*CF2*CKLNO (IPW1) 
CEPHN=CEPHN+DAH*DJP*CF2 *CI1*CKLNO (IPW1) 
ELSE 
CETHN=CETHN+DAH* DJP* CF2*CKLNO (IPW1) 
CE PHN=CEPHN+DAH*DJP*CF2*CI1*CKLNO (IPW1) 
END IF 
END IF 
TPOCRTHE .EO. PI) THEN 
iF {Neel 07) THEN 
CETHN=CETHN-DAH* DJP* CF1*CKLNO (IPW1 ) 
CE PHN=CEPHN+DAH* DJP*CF1*CI1*CKLNO ( IPW1) 
ELSE 
CETHN=CETHN+DAH* DJP*CF1*CKLNO (IPW1) 
CEPHN-CEPAN+DAH-“DIP*CF1*Cil*CKLNO(IPW1) 
END LF 
END IF 
CONTINUE 
ELSE 
DO 2200 IP=0, IQMAX 
PP=1P +1 
IPW1=4*IP+1 
IPW3=IPW1+2 
JP=MOD (IP, 4) 
CEl-Cit* oP 
CF2Z-CiZ7 oF 


DJP=DJPS (INP) 
DE nee EO. 0.) > THEN 

Pee Ne. et ee) THEN 
CETHN=CETHN-DAH* DJP* CF2* (CKLN (IPW1) -CI1*CKLN(IPW3) ) 
CEPHN=CEPHN+DAH* DJP*CF2* (CI1*CKLN (IPW1) +CKLN (IPW3) ) 

ELSE 
CETHN=CETHN+DAH* DJP*CF2* (CKLN (IPW1)+CI1*CKLN (IPW3) ) 
CEPHN=CEPHN+DAH*DJP*CF2* (CI1*CKLN(IPW1) -CKLN (IPW3) ) 


OD 


ENDSLE 
END IF 
FP OtRITHE ~EO. PIL) THEN 
TES ON  -GTo20). THEN 
CETHN=CETHN-DAH* DJP*CF1* (CKLN (IPW1)+CI1*CKLN (IPW3) ) 
CEPHN=CEPHN+DAH=DIP*CF1* (ClI1*CKLN(IPW1) —CKLN (I PW3)) 
ELSE 
CETHN=CETHN+DAH* DIP*CrI~ (CAEN (IPW1) -CIl*CKLN (LPW3)} 
CEPHN=CEPHN+DAH*DJP*CF1* (CI1*CKLN(IPW1)+CKLN (IPW3) ) 


END IF 
END IF 
2200 CONTINUE 
END IF 
3000 RETURN 
END 
Cc 


SUBROUTINE KLCRNT (DN, CKLN, CIW) 
CEE KKEKKKKKKKKEKEKKEKKKKKKKKKKKKEKKKKKKKKKKEKKKKKKKKKKKKKKKKKKKKEKKKKKKKKKKKKKKKK 
C This subroutine computes equivalent currents K and L on inner and outer 
C surfaces respectively. 
INCLUDE “REALTP. INC ' 
INCLUDE '‘CMPXTP.INC' 
INCLUDE ‘MAINDM.INC' 


(& 
DIMENSION CKLN(KCRNT) , CIW(KCRNT) 
DIMENSION CKLNP(KCRNT) , CKLNN (KCRNT) 
DIMENSION CRPQ31(KCRNT, KCRNT) , CRPQ32 (KCRNT, KCRNT) 
DIMENSION CKPHE (61.61) CKZ2P (6, 61),CLPHP(61, 61), 
= CLZP (Gi oy 4 CKPHN( OL, OL) -CKZN(GL, 62), 
+ CLPHN (617564), CLZN(61,61) 
COMMON /CKLMTX/ CKPHP,CKZP,CLPHP,CLZP, CKPHN, CKZN, CLPHN, CLZN 
COMMON /XPOTMP2/ CRPOQ31,CRPQ32 
COMMON /CRNTDM/ NMAX,MXNG, IOMAX, IOMAX1, IOMAX2, IXCRNT, ICRNT, MXQEG, 
+ MXOQOG : 
COMMON /INPUT3/ RTHE, RDELT, RPHI, RDELP,NTHTAO,NPHI, THESIN, THECOS, 
+ RHPI 
COMMON /INPUT4/ IZ,IK,IS,NYSM 
G 


DO 400 IX=1,ICRNT 
CKLNP (IX) =CZERO 
CKLNN (IX) =CZERO 

400 CONTINUE 
PHI1=DN*RPHI 
CPHI=CONE* COS (PHI1) +CLi*SIN(PHIL) 


eC 
DO 600 IX=1,ICRNT 
DOws00 TY=atCRNT 
ChmNP (i) =CKLNP (1X) +CRPO3S1(IX, LY) *CKEN( IY) 
500 CONTINUE 
600 CONTINUE 
c 
DO 800 IX=1,ICRNT 
DO 700 IY=1,ICRNT 
CKLNN (IX) =CKLNN (1X) +CRPQ32 (IX, IY) *CKLN (TY) 
700 CONTINUE 
800 CONTINUE 
Cc 


DOruGo0; Ir=-30,30 

DF=IF 

DPHI=DF*PI/32.D0 

DNPHI=DN* DPHI 

CPHI=CONE*DCOS (DNPHI ) +CI1*DSIN (DNPHT) 
Hom S00 J2=-30,30 

DIZ=J3Z 

DZ=DIZ/32.D0 

DV=DACOS (DZ) 

SINV=DSIN (DV) 
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CKPHP( IF, JZ) =CZERO 
CKZP(IF,JZ) =CZERO 
CLPHP (IF,JZ) =CZERO 
CLZP(IF,JZ) =CZERO 
CKPHN (IF,JZ) =CZERO 
CKZN (IF, JZ) =CZERO 
CLPHN(IF,JZ) =CZERO 
CLZN (IF,JZ) =CZERO 


DO 1400 IP=0,IQMAX 
P=IP 
P1=IP+1 
DPV=P* DV 
DP1V=P1*DV 
DCSPV=DCOS (DPV) 
SINP1V=DSIN(DP1V) 
IPX1=IP*4+1 
IPX2=IPX1+1 
IPX3=IPX2+1 
IPX4=IPX3+1 
CKPHP (IF, JZ) =CKPHP(IF,JZ)+CPHI* (CKLNP (IPX1)+CIW(IPX4) ) 

+ *DCSPV/PI/SINV 
CKZP(IF,JZ)=CKZP(IF,JZ)+CPHI* (CKLNP(IPX2) -CIW(IPX3))*SINP1V/PI 
CLPHP (IF,JZ)=CLPHP(IF,JZ)+CPHI* (CKLNP (IPX3) -CIW(IPX2) ) 

fe *DCSPV/PI/SINV 
CLZP(IF,JZ)=CLZP(IF,JZ)+CPHI* (CKLNP (IPX4) +CIW(IPX1)) *SINP1V/PI 
CKPHN (IF,JZ) =CKPHN (IF, JZ) +CPHI* (CKLNN (IPX1) -CIW(IPX4) ) 

+ *DCSPV/PI/SINV 
CKZN (IF,JZ)=CKZN(IF,JZ)+CPHI* (CKLNN (IPX2)+CIW(IPX3))*SINP1V/PI 
CLPHN(IF,JZ)=CLPHN(IF,JZ)+CPHI* (CKLNN(IPX3)+CIW(IPX2) ) 

= *DCSPV/PI/SINV 
CLZN (IF,JZ) =CLZN(IF,JZ) +CPHI* (CKLNN(IPX4) -CIW(IPX1)) *SINP1V/PI 

1400 CONTINUE 


sy 9) CONTINUE 
1600 CONTINUE 
RETURN 
END 


J. SUBROUTINE RCSPAREA 


SUBROUTINE RCSPAREA(CESTH, CESPH) 
CER KEK RRR KEK KKK KKK KKK KK KR KK KEK RK KK KKH KKK RK KKK KEK RK KEK KKK KKK KK KEK KK KKK KKK KKK KK 
C This subroutine computes cross section per projected area in all direction. 
c 

INCLUDE ‘REALTP.INC' 

INCLUDE 'CMPXTP.INC' 


S 
REAL*4 ARCSPPA1,ARCSPPA2,APHASE1, APHASE2 
COMMON /GCONST/ DKH, DKA, HH,HA,HSQ, DASQ, HSQN, DASON, HHSQ, HDASQ, 
a RAH, RAHSQ, DHA, DAH 
COMMON /INPUT3/ RTHE, RDELT, RPHI, RDELP,NTHTAO, NEFH1, THESIN, THECOS, 
Da RHPI 
c 


ESTH=ABS (CESTH) 
ESPH=ABS (CESPH) 
ESTHSQO=ESTH* ESTH 
ES PHSQ=ESPH*ESPH 
TR (RTHE EO] REPL) THEN 
RCSPPA1=PI* ESTHSQ/DKH/ DKA 
RCSPPA2=PI* ESPHSQ/DKH/ DKA 
BESE 1h ((RTHE .E0O. ZERO) .OR. (RTHE EO. Pi) Teen 
RCSPPA1=4.D0*ESTHSQ/DASQ 
RCSPPA2=4 .D0* ESPHSQ/DASQ 
ELSE 


ul 


AP=ABS (QUAR* DASQ* THECOS ) +ABS (DKH* DKA* THESIN/ PT) 
RCSPPA1=ESTHSQ/AP 
RCSPPA2=ESPHSQ/AP 
END IF 
ARCSPPA1=RCSPPA1 
ARCSPPA2=RCSPPA2 


ER1=REAL (CESTH) 

EIT1=IMAG (CESTH) 

ER2=REAL (CESPH) 

EI2=IMAG (CESPH) 

PHASE1=DATANZ2 (EI1, ER1) 

PHASE2=DATANZ2 (EI2, ER2) 

APHASE1=PHASE1 

APHASE2=PHASE2 

WRITE (21,*) ARCSPPA1,APHASE1,ARCSPPA2,APHASE2 


RETURN 
END 


AE 
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